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Introduction 



QCD at Hadron Colliders 

Quantum Chromodynamics is the accepted theory of the strong interactions. After 
decades spent to test this theory and after its successfull description of the rich and 
complex phenomenology of the hadronic world, nobody seriously doubts that the simple 
Lagrangian which is at the core of its formulation is the correct one. 

At the same time, the wider dynamics of what has come to be known as the Standard 
Model of elementary particles has been tested with an incredible accuracy at LEP and 
at the Tevatron, confirming the validity of the formulation of the model as a local gauge 
theory based on the gauge group SU (3) x SU (2) x U (l)y , where U (l)y is the abelian group 
describing the hypercharge. After the incorporation of spontaneous symmetry breaking of 
the SU{2) x i7(l) symmetry to U{l) em via the Higgs mechanism, the model can describe 
the electroweak interactions, leaving a broken gauge theory which is exact only in the 
color sector (S77(3)) and in the electromagnetic U(l) em abelian group. The application of 
these field theoretical methods to the analysis of specific processes is a difficult scientific 
enterprise and requires the identification of special theoretical tools, the proof of theorems 
in quantum field theory and the development of sophisticated software - based on theory 
- for accurate predictions of the experimental observables of a large array of processes. 
In this effort one discovers that there are energy regions of the interactions where a 
perturbative picture of the theory can be slowly built, thanks to the property of asymptotic 
freedom, and others where this is not possible due to a strong coupling constant. In 
the perturbative region, being the theory still characterized by the property of quark 
confinement, a class of theorems known as factorization theorems allow the definition of a 
calculable scheme within which to compare the theory with the experiment. The essential 
ingredient of this procedure is the possibility to separate - in inclusive processes - the 
perturbative part of a process - what is known as the hard scattering - from the non- 
perturbative parts, termed parton distribution functions. In exclusive processes, at large 
momentum transfers, a similar picture holds, and allows to describe nucleon form factors 
and other elastic processes by hadronic wave functions, which nowadays are special cases of 
some non-local matrix elements on the light cone termed nonforward parton distributions. 
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As we have just mentioned, perturbative QCD appears naturally at large energy and 
momentum transfers and is light-cone dominated. Parton distributions are non-local 
correlators defined for light cone separations of the field operators and their definition in- 
trinsically introduces a scale. This description goes under the name of parton model. The 
changes induced on this operators due to the changes of the scale at which they are de- 
fined can be described by renormalization group methods. In the case of ordinary parton 
distributions these equations are called DGLAP equations and they have been studied in 
various contexts (for polarized and unpolarized collision processes) and at various orders 
in perturbation theory. One of the objective of this thesis is to elaborate on the phe- 
nomenology of these equations in great detail, to next-to-leading order (NLO) (chapter 
1) and up to next-to-next-to-leading order (NNLO) in a s , the QCD strong coupling con- 
stant (chapter 2). This perturbative order is likely to be the "final state" for the study of 
the DGLAP equations in perturbative QCD, since it provides an accuracy which clearly 
satisfies the requirements of precisions needed at hadron colliders. The study of these 
equations is rather comprehensive, involves all the leading-twist distributions up to NLO 
and is further extended to NNLO for the unpolarized distributions. Applications of these 
studies include a prediction for the total cross section for Higgs detection at the LHC at 
NNLO (chapter 3), where the dependence on the predictions on the factorization scale is 
studied in depth and an NLO application for the study of the Drell Yan process for the 
PAX experiment (chapter 1) at few GeV's. These applications will all be discussed in the 
first three chapters of this thesis work. In chapter 4 and in chapter 5 we elaborate on a 
kinetic approach to the DGLAP evolution by introducing a Kramers- Moyal expansion of 
this equation after testing both analytically and numerically the positivity of the kernels 
up to NLO. In chapter 5 we continue and enlarge this analysis and apply it to the study 
of an inequality due to Soffer and to its behaviour under the action of the renormalization 
group. Incidentally, we also present an ansatz for the evolution of nonforward parton 
distributions and prove the existence of a kinetic form in the nonsinglet case. 

QCD in the Astroparticle Domain 

In the second part of this thesis we will try to provide some applications of standard 
perturbative QCD methods in the study of the extensive air showers which are generated 
by proton and neutrino primaries in high and ultra high energy cosmic rays (UHECR). 
QCD enters in the modeling of the fragmentation region and our QCD predictions are 
appropriately interfaced with an air shower simulator which is well known in the cosmic 
rays community, CORSIKA. 

This second part of the thesis can be read quite independently from the first part, but 
has been made possible by the developments of various computational tools described in 
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the first part. In this second part we will define some observables of the air showers which 
can help us discriminate between standard scenarios and "new physics" scenarios in the 
context of UHECR (chapter 6), having in mind the ongoing experimental efforts of various 
experimental collaborations, such as Auger, to assess the existence (or the absence) of a 
cutoff in the upper part of the cosmic rays spectrum. These issues have been addressed 
in two papers. In a first paper we have introduced some specific observables for the 
study of the extensive air showers and modified the hadronization codes used by the 
various CORSIKA routines and studied their impact on the lateral distributions and the 
multiplicity distributions, while in a second paper, on which chapter 7 is based, we have 
analyzed in great detail the characteristic patterns of the showers in the context of one 
specific scenarios for new physics which goes under the name of theory of extra dimensions. 
In this last case an attempt has been made to understand the so called Centauro events 
in cosmic rays physics as evaporating mini black holes, which are predicted to form in 
theories with extra dimensions. The methodology that we develope is quite general and 
can be used to discern between standard and new physics scenario also in the context 
of other theories, since we propose and exemplify a general strategy. This second part 
of the thesis is heavily computational and the author has developed software for the 
complex analysis of the data generated by the simulator. In particular, the modeling of 
the hadronization of the evaporating mini black hole has been developed using a basic 
physics picture of the stages appearing in its decay and a hadronization pattern of the 
partonic states which is in concordance with the "democratic" coupling of gravity to all 
the states of the standard model. 
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Chapter 1 



Direct solution of Renormalization 
Group Equations of QCD in x-space: 
NLO implementations at leading twist 

1.1 Objectives of the chapter and Background 

In this chapter we formulate and implement an algorithm for the solution of evolution 
equations in QCD which has some specific peculiarities compared to other methods based 
on the use of Mellin moments or of Laguerre polynomials. We present the motivation of 
this work and discuss the implementation of the method up to next to leading order (NLO) 
in a s , the QCD strong coupling constant. We also briefly discuss some of the technical 
aspects of the derivation of the recursion relations which allow to rewrite an operatorial 
solution of these equations in terms of some special functions determined recursively. The 
method is implemented numerically and a brief discussion of the structure of the program 
is included. The content of this chapter is based on the original research article [25]. In 
the final part of this chapter we provide an application of the results on the evolution of 
the transverse spin distributions to the case of Drell-Yan lepton pair production in the 
few GeV region, near the J/ip resonance, up to NLO in a s , where an enhancement is 
expected. Predictions for the transverse spin asymmetries of this process are presented. 

1.1.1 Factorization and Evolution 

Factorization theorems play a crucial role in the application of perturbation theory to 
hadronic reactions. The proof of these theorems and the actual implementation of their 
implications has spanned a long time and has allowed to put the parton model under a 
stringent experimental test. Prior to embark on the discussion of our contributions to the 
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1.1. Objectives of the chapter and Background 



study of evolution algorithms in Bjorken x-space, we provide here a brief background on 
the topic in order to make our treatment self-contained. 

In sufficiently inclusive cross sections, leading power factorization theorems allow to 
write down a hadronic cross section in terms of parton distributions and of some hard 
scatterings, the latter being calculable at a given order in perturbation theory using the 
fundamental QCD Lagrangian. Specifically, for a hadronic cross section, for instance a 
proton-proton cross section a pp , the result of the calculation can be summarized by the 
formula 

a PP = J2 f dx i I dx2fh 1 ^f(xi,Q 2 )f h2 ^ f (x 2 ,Q 2 )a(x 1 ,X2,s,i,Q 2 ) , (1.1) 
J Jo Jo 

where the integral is a function of some variables X\ and x 2 which describe the QCD 
dynamics at parton level in the Deep Inelastic Scattering (DIS) limit or, equivalently, at 
large energy and momentum transfers. These variables are termed Bjorken variables and 
are scale invariant. This formula is a statement about the computability of a hadronic 
collisions in terms of some "building blocks" of easier definition. 

The variable Q 2 , in the equation above, can be identified with the factorization scale 
of the process, a can be computed at a given order in perturbation theory in an expansion 
in a s , while the fh t ^f(x,Q 2 ) are the parton distributions. These describe the probability 
for a hadron h to prepare for the scattering a parton /, which undergoes the collision. 
An equivalent interpretation of the functions fhi^f(x,Q 2 ) is to characterize the density 
of partons of type / into a hadron of type h. A familiar notation, which simplifies the 
previous notations shown above, is to denote by qi(x, Q 2 ) the density of quarks in a hadron 
(a proton in this case) of flavour i and by g(x, Q 2 ) the corresponding density of gluons. 
For instance, the annihilation channel of a quark with an antiquark in a generic process 
is accounted for by the contribution 

J dxi dx2q(x 1 ,Q 2 )q(x 2 ,Q 2 )(t q q(x 1 ,X2,Q 2 ) (1.2) 

and so on for the other contributions, such as the quark-gluon sector (qg) or the gluon- 
gluon sector (gg) each of them characterized by a specific hard scattering cross section 
a qg , or a qg , and so on. In this separation of the cross section into contributions of hard 
scatterings a and parton distributions f(x,Q 2 ) the scale at which the separation occurs 
is, in a way artificial, in the sense that the hadronic cross section a should not depend on 
Q 2 or 

da 

5F = °- (L3) 

However, a perturbative computation performed by using the factorization formula, how- 
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ever, shows that this is not the case, since the perturbative expansion of a 

a = a^ + a s (Q 2 ) 2 aW+a s (Q 2 )aV (1.4) 

naturally brings in the dependence on the factorization scale Q. This dependence is weaker 
if we are able to push our computation of the hard scattering to a sufficiently high order 
in a s . The order at which the perturbative expansion stops is also classified as a "leading 
order" (LO), "next-to-leading order" (NLO) or, even better, a "next-to-next-to-leading" 
(NNLO) contribution if more and more terms in the expansion are included, 

At the same time, the parton distributions f(x, Q 2 ) carry a similar dependence on the 
scale Q, which is summarized in a renormalization group equation (RGE) which resums 
the logarithmic violations to the scaling behaviour induced by the perturbative expansion. 
Also in this case we need to quantify this effect and reduce its impact on the prediction of 
the cross section. The topic of this chapter, therefore, is about the quantification of this 
effect and the implementation of an algorithm which reorganizes the solutions of these 
RGE's in a suitable way to render the numerical implementation quick and accurate up 
to NLO. This strategy is fully worked out for all the leading twist distributions, with 
and without polarization. The results of these chapter can be applied for the precise 
determination of both polarized and polarized cross sections in pp collisions. In chapter 
3 this methodology will be pushed even further, to NNLO, which is the state of the art 
in QCD and will be applied to the case of the Higgs croiss section, which is important for 
the discovery of this particle at the LHC. 

1.2 Parton Dynamics at NLO: a short overview 

Our understanding of the QCD dynamics has improved steadily along the years. In fact we 
can claim that precision studies of the QCD background in a variety of energy ranges, from 
intermediate to high energy - whenever perturbation theory and factorization theorems 
hold - are now possible and the level of accuracy reached in this area is due both to 
theoretical improvements and to the flexible numerical implementations of many of the 
algorithms developed at theory level. 

Beside their importance in the determination of the QCD background in the search of 
new physics, these studies convey an understanding of the structure of the nucleon from 
first principles, an effort which is engaging both theoretically and experimentally. 

It should be clear, however, that perturbative methods are still bound to approxima- 
tions, from the order of the perturbative expansion to the phenomenological description 
of the parton distribution functions, being them related to a parton model view of the 
QCD dynamics. 
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1.2. Parton Dynamics at NLO: a short overview 



Within the framework of the parton model, evolution equations of DGLAP-type - and 
the corresponding initial conditions on the parton distributions - are among the most im- 
portant blocks which characterize the description of the quark-gluon interaction and, as 
such, deserve continuous attention. Other parts of this description require the computa- 
tion of hard scatterings with high accuracy and an understanding of the fragmentation 
region as well. The huge success of the parton model justifies all the effort. 

In this consolidated program, we believe that any attempt to study the renormalization 
group evolution describing the perturbative change of the distributions with energy, from 
a different - in this case, numerical - standpoint is of interest. 

In this chapter we illustrate an algorithm based on the use of recursion relations for 
the solution of evolution equations of DGLAP type. We are going to discuss some of the 
salient features of the method and illustrate its detailed implementation as a computer 
code up to next-to-leading order in a s . 

In this context, it is worth to recall that the most common method implemented so 
far in the solution of the DGLAP equations is the one based on the Mellin moments, with 
all its good points and limitations. 

The reason for such limitations are that, while it is rather straightforward to solve 
the equations in the space of moments, their correct inversion is harder to perform, since 
this requires the computation of an integral in the complex plane and the search for an 
optimized path. 

In this respect, several alternative implementations of the NLO evolution are available 
from the previous literature, either based on the use of "brute force" algorithms [80] or on 
the Laguerre expansion [63, 45], all with their positive features and their limitations. 

Here we document an implementation to NLO of a method based on an ansatz [109] 
which allows to rewrite the evolution equations as a set of recursion relations for some scale 
invariant functions, A n (x) and B n (x), which appear in the expansion. The advantage, 
compared to others, of using these recursion relations is that just few iterates of these are 
necessary in order to obtain a stable solution. The number of iterates is determined at run 
time. We also mention that our implementation can be extended to more complex cases, 
including the cases of nonforward parton distributions and of supersymmetry. Here we 
have implemented the recursion method in all the important cases of initial state scaling 
violations connected to the evolutions of both polarized and unpolarized parton densities, 
including the less known transverse spin distributions. 

Part of this chapter is supposed to illustrate the algorithm, having worked out in some 
detail the structure of the recursion relations rather explicitly, especially in the case of 
nonsinglet evolutions, such as for transverse spin distributions. 

One of the advantages of the method is its analytical base, since the recursion relations 
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can be written down in explicit form and at the same time one can perform a simple 
analytical matching between various regions in the evolutions, which is a good feature of 
x-spaced methods. While this last aspect is not relevant for ordinary QCD, it is relevant 
in other theories, such as for supersymmetric extensions of the parton model, once one 
assumes that, as the evolution scale Q raises, new anomalous dimensions are needed in 
order to describe the mixing between ordinary and supersymmetric partons. 



1.3 Definitions and Conventions 

In this section we briefly outline our definitions and conventions. 

We will be using the running of the coupling constant up to two-loop level 



a s (Q 2 ) ' 7r 



A> log(Q 2 /A|^) 



1 _ A iogiog(Q7Ak) 4 () 



(1.5) 



where 



and 



A> = y Nc ~ t Tf > Pi = Y N c-Y Ncn '- 2CFn f> (L6) 

N 2 — 1 4 1 
^c = 3, C F = ^- = -, T f = T R n f = -n f , (1.7) 

where Nc is the number of colors, n/ is the number of active flavors, which is fixed by the 
number of quarks with m q < Q. We have taken for the quark masses m c — 1.5 GeV, mj = 
4.5 GeV and m t = 175 GeV, these are necessary in order to identify the thresholds at 
which the number of flavours rif is raised as we increase the final evolution scale. 
In our conventions Aqcd is denoted by A^-? and is given by 

A |M,5,6) = a24g) o 200, 0.131, 0.050 GeV. (1.8) 

We also define the distribution of a given helicity (±), / ± (a;, Q 2 ), which is the prob- 
ability of finding a parton of type / at a scale Q, where / = qi,qi,g, in a longitudinally 
polarized proton with the spin aligned (+) or anti-aligned (— ) respect to the proton spin 
and carrying a fraction x of the proton's momentum. 

As usual, we introduce the longitudinally polarized parton distribution of the proton 

Af(x,Q 2 ) = f + (x,Q 2 )-f-(x,Q 2 ). (1.9) 



We also introduce another type of parton density, termed transverse spin distribution, 
which is defined as the probability of finding a parton of type / in a transversely polarized 
proton with its spin parallel (|) minus the probability of finding it antiparallel (J,) to the 
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proton spin 

A T f(x, Q 2 ) = f\x, Q 2 ) - fL(x, Q 2 ). (1.10) 
Similarly, the unpolarized (spin averaged) parton distribution of the proton is given 

by 

f(x, Q 2 ) = f+(x, Q 2 ) + f-(x, Q 2 ) = f(x, Q 2 ) + f l (x, Q 2 ). (1.11) 

We also recall that taking linear combinations of Equations (1.11) and (1.9), one recovers 
the parton distributions of a given helicity 

sHx ^ )= nam^m. (1 . 12) 

In regard to the kernels, the notations P, AP, A^P, P ± , will be used to denote the 
Altarelli-Parisi kernels in the unpolarized, longitudinally polarized, transversely polarized, 
and the positive (negative) helicity cases respectively. 

The DGLAP equation is an integro-differential equation whose general mathematical 
structure is 

^-^f(x, Q 2 ) = ®s(Q 2 )) ® f(x, Q 2 ), (1.13) 
where the convolution product is defined by 

[a ® b] (x) = f ^-a (*) b(y) = f %(y)b (*) . (1.14) 

Jx y \yj Jx y \y J 
Let us now turn to the evolution equations, starting from the unpolarized case. Defin- 



ing 

# } = ft ± ft , = E # } , X, = # } " -<7 {+) , 

n f 



(1.15) 



the evolution equations are 

d 



dlogQ 2 
d 



qt\x, Q 2 ) = P NS - (x, a s (Q 2 )) ® qt\x, Q 2 ), (1.16) 
^(x,Q 2 ) =P N s+(x,® s (Q 2 ))®k(x,Q 2 ), (1-17) 



dio g g 

for the nonsinglet sector and 

d / q (+) (x,Q 2 ) \ = ( P qq (x,a s (Q 2 )) P qg (x,a s (Q 2 )) \ ( q^(x,Q 2 ) 
dlogQ 2 ^ g(x,Q 2 ) ) \P gq (x,a s (Q 2 )) P gg (x,a s (Q 2 )) ) { g(x,Q 2 ) 

(1.18) 

for the singlet sector. 
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Equations analogous to (1.15 - 1.18), with just a change of notation, are valid in the 
longitudinally polarized case and, due to the linearity of the evolution equations, also for 
the distributions in the helicity basis. In the transverse case instead, there is no coupling 
between gluons and quarks, so the singlet sector (1.18) is missing. In this case we will 
solve just the nonsinglet equations 

d A T qt\x,Q 2 ) = A T P NS -(x,a s (Q 2 )) <g> A T qt\x,Q 2 ), (1.19) 



2 



dlogQ 

^±^A T qt\x,Q 2 ) = A T P NS+ (x,a s (Q 2 ))^A T qi + \x,Q 2 ). (1.20) 

We also recall that the perturbative expansion, up to next-to-leading order, of the 
kernels is 

P(x, a.) = (g) P<°> (x) + 2 PW (x) + . . . . (1.21) 

1.4 Mathematical structure of the kernel 

Altarelli-Parisi kernels are defined as distributions. The most general form is the following 

P{x) = p x {x) + -^L + p 3 S(i - X ), (1.22) 
I 1 - x ) + 

with a regular part Pi (x) , a "plus distribution" part P 2 (x) / ( 1 — x) + and a delta distribution 
part P3<5(1 — x). For a generic function a(x) defined in the [0, 1) interval and singular in 
x — 1, the plus distribution [a(x)]+ is defined by 

f 1 f(x)[a(x)] + dx = f 1 (/Or) - /(l)) a(x)dx, (1.23) 

where /(x) is a regular test function. Alternatively, an operative definition (that assumes 
full mathematical meaning only when integrated) is the following 

[a(x)]+ = a(x) - 5(1 - x) C a(y)dy. (1.24) 

Jo 

From (1.24) it follows immediately that each plus distribution integrate to zero in the 
[0, 1] interval 

/ [a(a:)] + da: = 0. (1.25) 
Jo 



Convolution of a generic kernel. 
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1.4. Mathematical structure of the kernel 



We want to make the convolution of the generic kernel (1-22) with a function f(x). 
To improve numerical stability, in the code we indeed multiply the functions to be 
convolved with the kernel (the coefficients A n and B n ) by an x factor. To compute 
the formulas that we implement in the code, we introduce the notation f(x) = xf(x). 
The treatment of the regular and the delta-function parts is trivial 

P^x) ® f(x) = xP x {x) ® f(x) = x f 1 ^Pi(y)/ (-) = C &yPi(y)f (-) (1.26) 

Jx y \yj Jx \yj 

P 3 S(1 ~ x) ® f(x) = ^ ^-P 3 5(l - y)f (-) . (1.27) 
Jx y \yj 

Let us now treat the more involved case of the plus distribution part 

Pz(x) £l s P 2 (x) „ , f f 1 dy 



Jx y l-y \y ) Jo 1 - y Jx y 



x 



Jx y l-y \y Jo 1 - 



dy 

V \V J ■ Jo l-y 

Jx y l-y \yj Jx 



i-y 



which yields 



(1-x). 



-P 2 (l)f(x) £ ^- (1.28) 

/^^^^-^^^V/jMl-,), (1.29) 
Jx y 1 — y 



m = /' + /» io g(1 - (i.30) 



Mellin Moments 

The n-th Mellin moment of a function of the Bjorken variable f(x) is defined by 



/„= [ 1 x n - 1 f(x)dx. (1.31) 

J 
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An important property of Mellin moments is that the Mellin moment of the convolu- 
tion of two functions is equal to the product of the individual Mellin moments 

[f®g] n = fng n . (1-32) 

Let us prove it. 

[/«$]„ = f dxx n ~ l [f®g](x) = f'dxx"- 1 [ 1( ^f(y)g(*). (1.33) 
Jo Jo Jx y \y J 

Exchanging the x and y integrations 

[/««]„ = (1-34) 
and introducing the new variable z = x/y 

[f®g] n = j o 1 dyf(y)y n ^j Q 1 dzz n - 1 g(z) = f n g n . (1.35) 

This leads to an alternative formulation of DGLAP equation, that is also the most 
widely used to numerically solve the evolution equations. By taking the first Mellin 
moment of both sides of the integro-differential equation (1.13) we are left with the 
differential equation 

Z] ^/i(Q 2 )=Pi(Q 2 )/i(Q 2 ) (1-36) 

that can be easily solved to give 

fi(Q 2 )= C f(x,Q 2 )dx. (1.37) 
Jo 

To get the desired solution f(x,Q 2 ) there is a last step, the inverse Mellin transform 
of the first moment of the parton distributions, involving a numerical integration 
on the complex plane. This is the most difficult (and time-consuming) task that the 
algorithms of solution of DGLAP equation based on Mellin transformation - by far the 
most widely used - must accomplish. But the most severe restriction of the flexibility 
of such evolution codes (for example QCD-Pegasus by Vogt [116]) as compared to 
x-space methods is that the initial distribution are needed in a form which facilitates 
the inversion of the moments, i.e. a functional form; but very often parton sets are 
given in discrete grids of x and Q values. 



From Feynman diagrams calculations one can get just the regular part Pi(x) of each 
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kernel. The remaining distributional parts (plus distribution and delta distribution) 
emerge from a procedure of regularization, that introduce the plus distribution part to 
regularize the eventual singularity in x — 1 and the delta distribution to fulfill some 
physical constraints, the sum rules. 

The first one is the baryon number sum rule (BNSR), asserting that the baryon number 
(number of quarks less number of antiquarks) of the hadron must remain equal to its initial 
value (3 in the case of the proton) throughout the evolution, i.e. for each value of Q 2 



qt\Q 2 )= f 1 q { -\x,Q 2 )dx = 3. 
Jo 



(1.38) 



Deriving (1.38) with respect to logQ 2 and having in mind that q( ) evolves with P^ s , we 
get 

(1.39) 



dx 



(x) = 0. 



PlsiQ 2 )®<i { -\Q 2 )_ 

Making use of the property of the Mellin moment of a convolution (1.32) this implies 

QT 1 P% s (x,Q 2 )dx) QTV-JfoQ^dx) = 0, (1.40) 

from which, using (1.38), we find the BNSR condition on the kernel 

/' P% s (x)dx = 0. (1.41) 
Jo 

The other constraint is the momentum sum rule (MSR), asserting that the total mo- 
mentum of the hadron is constant throughout the evolution. Having in mind that x is 
the fraction of momentum carried out by each parton, this concept is translated by the 
relation 

J (xq {+) (x, Q 2 ) + xg(x, Q 2 )) dx = 1 (1.42) 

that must hold for each value of Q 2 . Deriving with respect to logQ 2 and using the singlet 
DGLAP equation 



+ 



P q9 {Q 2 ) ® g{Q 2 )] (x) 
P gg {Q 2 )®g{Q 2 )]{x)} = o. 



I' dxx{[P qq (Q 2 )®qt + \Q 2 )] (x) 

J 

+ [p 9q {Q 2 )®q {+ \Q 2 )_ 

Using (1.32) we get 

£ x (P qq (x, Q 2 ) + P gq (x, Q 2 )) dx\ [jf 1 xq^ (x, Q 2 )dx 
£ x (P qg (x, Q 2 ) + P gg (x, Q 2 )) dx jT 1 xg(x, Q 2 )dx 



(1.43) 



+ 



0, 



(1.44) 
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from which we find the MSR conditions on the singlet kernels 

£ x [P qq (x, Q 2 ) + P gq (x, Q 2 )) dx = 0, 

x (P qg (x,Q 2 ) + P 99 (x,Q 2 )) dx = 0. 



(1.45) 
(1.46) 



The kernels regularization procedure: an example. 

We illustrate now an example of the regularization procedure of the Altarelli-Parisi 
kernels through the sum rules. The LO kernels computed by diagrammatic techniques 
for x < 1 are 



P£\x) = P$l{x)=C F 



l + x 2 



1 — x 



1 — x 



-l-x 



P^(x) = 2T f 
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P£\x)=C F 



x 2 + (l-x) 2 



l + (l-x) 



X 



Pj®(x)=2N c 



1 +--2 + x(l-x) 



(1.47) 

(1.48) 
(1.49) 

(1.50) 



.1 — x x 

We want to analytically continue this kernels to x — 1 curing the ultraviolet singular- 
ities in (x) and P^ (x) . We start introducing the plus distribution prescription in 
PqqK x )- We make the replacement 



(1.51) 



l-x 



(l-x) 



+ 



to avoid the singularity and we add a term k8(l — x) (where k has to be determined) 
to fulfill the BNSR (1.41). So we have 



(l-x). 



- 1 - x + kS(l-x) 



(1.52) 



Imposing by the BNSR that P£\x) integrates to zero in [0, 1] and remembering that 
the plus distribution integrates to zero we get 



Jo P ^ ix)d 



x = C, 



0, 



(1.53) 
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hence k = 3/2, and the regularized form of the kernel is 



(l-x) + 



l-x + -5(1 -x) 



(1.54) 



Noticing that 



f 1 x n f 1 x - 1 + 1 , f 1 
/ 7 ^ax = / — ax — . 

Jo (1 — x) + Jo (1 — x)+ Jo \ (1 — x). 



-1 + 



dx = -1 (1.55) 



it can be easily proved that the MSR (1.45) is satisfied. Let us now regularize P gg (x). 
We make the replacement 



(1 — x)+ X 
Imposing the other MSR (1.46) we get 



11 / 

H 2 + x(l-x) 



+ kS(l-x). (1.56) 



2N, 



a 



from which we find 



x 



:i-x). 



+ 1 -2x + x 2 (l -x) 



+ kx5(l — x) 



+27) [x 3 +x(l -x) 2 ]}dx = 0, (1.57) 



k = T N °-l T > 



Po 
2 ' 



so the regularized form of the kernel is 



PM{X) = 2N C 



11 / 

- + --2 + X (1-x) 

(1 — X) + X 



+ !*(!-*). 



(1.58) 



(1.59) 



1.5 The Ansatz and some Examples 

In order to solve the evolution equations directly in x-space, we assume solutions of the 
form 

f(x, q 2 ) = £ ^ log" + a s(Q 2 ) £ log" ^4S ? (1-60) 



n=0 



Us{Ql) 



n=0 



for each parton distribution /, where Q defines the initial evolution scale. The justifi- 
cation of this ansatz can be found, at least in the case of the photon structure function, 
in the original work of Rossi [109], and its connection to the ordinary solutions of the 
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DGLAP equations is most easily worked out by taking moments of the scale invariant 
coefficient functions A n and B n and comparing them to the corresponding moments of 
the parton distributions, as we are going to illustrate in section 1.7. The link between 
Rossi's expansion and the solution of the evolution equations (which are ordinary differ- 
ential equations) in the space of the moments up to NLO will be discussed in that section, 
from which it will be clear that Rossi's ansatz involves a resummation of the ordinary 
Mellin moments of the parton distributions. 

Setting Q = Q in (1.60) we get 

f(x, Ql) = A (x) + a s (Q 2 )B (x). (1.61) 

Inserting (1.60) in the evolution equations, we obtain the following recursion relations for 
the coefficients A n and B n 

An +1 (x) = -^P (0 \x) A n (x), (1.62) 
Po 



B n+1 (x) = -B n (x) - J±-A n+1 (x) - ^P (0 \x) ® B n (x) - \p {1 \x) ® A n (x) (1.63) 

4ttPo Po ^Po 

obtained by equating left-hand sides and right-hand-side of the equation of the same 
logarithmic power in log™ a s (Q 2 ) and a s log™ a s (Q 2 ). Any boundary condition satisfying 
(1.61) can be chosen at the lowest scale Qq and in our case we choose 

B o (x) = 0, f(x,Ql)=A (x). (1.64) 

The actual implementation of the recursion relations is the main effort in the actual 
writing of the code. Obviously, this requires particular care in the handling of the singu- 
larities in x-space, being all the kernels defined as distributions. Since the distributions 
are integrated, there are various ways to render the integrals finite, as discussed in the pre- 
vious literature on the method [47] in the case of the photon structure function. In these 
previous studies the edge-point contributions - i.e. the terms which multiply 5(1 — x) in 
the kernels - are approximated using a sequence of functions converging to the 5 function 
in a distributional sense. 

This technique is not very efficient. We think that the best way to proceed is to 
actually perform the integrals explicitly in the recursion relations and let the subtracting 
terms appear under the same integral together with the bulk contributions (x < 1) (see 
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also [75]). This procedure is best exemplified by the integral relation 



dy 



x 



* y 



1 dy f(x/y) - yf(x) 



x y 



+ /(a;)log(l-x) 



(1.65) 



in which, on the right hand side, regularity of both the first and the second term is explicit. 
For instance, the singlet evolution equations become in the unpolarized case 



dgW(x) 
dlogQ 2 



= c. 



1 dy q( + \x/y) — yq^ + \x) 



x y 



+2T R n f 



+2qW(x) log(l - x) - ^ ^(1 + y)q™ (-) + \q 

Jx II \ II I z 

1 



(+), 



y 2 + (i-yY 



(1.66) 



dg(x) 

dio g g 2 



i dy 1 + (1 - 



+2iV c 



1 dyg(xfy) - y#(x) 



+ g(x) log(i - x) 



+ 



, Po / \ 



(1.67) 



1.6 An Example: The Evolution of the Transverse Spin 
Distributions 

LO and NLO recursion relations for the coefficients of the expansion can be worked out 
quite easily. We illustrate here in detail the implementation of a nonsinglet evolutions, 
such as those involving transverse spin distributions. For the first recursion relation (1.62) 
in this case we have 



A± +1 (x) = --A T pW(x)®A±(x) = 
Po 



Cf "a. 



1 dyyA±(y) - xA±(x) 



C F 



x y 
*dy 
x y 



y-x 
A$(y))+C F 



+ A±(x) log(l - x) 
2 \ 3 



+ 



Po 



(1.68) 



As we move to NLO, it is convenient to summarize the structure of the transverse kernel 
A T P^(x) as 
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A T P^\x) = K±(x)5(l -x) + Kf{x)S 2 {x) + Kf{x) \og{x) 



+Kf(x) log 2 (x) + Kf(x) \og{x) log(l -x) + Kf(x) 



Hence, for the (+) case we have 



+ K±(x). (1.69) 



A T P+.«(x) ® A+(x) = KfA+(x) + f *L [tf+(*)S 2 (z) + X 3 +(z) log(s) 

Jx y L 

+ l0g 2 (2)A-+( 2 ) + l0g(j) l0g(l - z)if 5 + (2)] At(y) + 

h7 1 jg g^M^M + A+n(x) log(1 _ J + /ff y (L70) 

y y — x j Jx y 



where z = x/y. For the (— ) case we get a similar expression. 
For the B^ +1 (x) we get (for the (+) case) 



b: +1 (x) 



-B:(x) + ^\2C F 



dyyA+(y)-xA+{x) 



x y 



y-x 



+ A+(x) log(l 



x) 



+ 



2 " I 4ttA) 
+ K+(z) \og{z) + log 2 (z)X 4 + (z) + log(z) log(l - z)tf 5 + (z) 



1 



4l(y) + 



dyyA+(y) -xA+(x) 



x y 



y-x 



+ A+(x) log(l-x) 



Jx y ) 



C f 



C r 



! o) 



- 1 dyyB±{y) - xB±(x) 



+ S±(x)log(l-x) 



-V/ 

where in the (+) case we have the expressions 



+ 



x 3 +(x) 

^ 4 +(*) 



—^(-2/1/(3 + 4tt 2 ) + iV c (51 + 44tt 2 - 216C(3)) + 9C F (3 - 4tt 2 + 48((3)) 

2C F (-2C F + N c )x 
1 + x 

C F {9C F -11N C + 2n f )x 



3(x- 1) 



C F NqX 

1-x 
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Kt{x) 



AC 2 F x 
1-x 



K U X ) = ~C F (10n f + N c (-Q7 + 3tt 2 )) 

y 

Kt(x) = l-C F (10n f + N c (-Q7 + 3tt 2 )), 

y 



(1.71) 



and for the (— ) case 

Ki(x) = ^C F (-2n f (3 + 4tt 2 ) + N c (51 + Un 2 - 216C(3)) + 9C F (3 - 4tt 2 + 48((3)) 
JZ _. 2C F (+2C F - N c )x 

K 2 {X) = 



K~{x) 
K7(x) 



l + x 

C F {9C F -llN c + 2n f )x 
3(x - 1) 

C F N c x 



Kz(x) = ~C F (10n f + N c (-Q7 + 3tt 2 )) 

Kj(x) = ~C F {10n f - 18C F (x - 1) + N c (-76 + 3tt 2 + 9x)). 

y 



(1.72) 



The terms containing similar distribution (such as "+" distributions and 5 functions) have 
been combined together in order to speed-up the computation of the recursion relations. 



1.7 Comparisons among Moments 

It is particularly instructive to illustrate here briefly the relation between the Mellin 
moments of the parton distributions, which evolve with standard ordinary differential 
equations, and those of the arbitrary coefficient A n (x) and B n {x) which characterize 
Rossi's expansion up to next-to-leading order. This relation, as we are going to show, 
involves a resummation of the ordinary moments of the parton distributions. 

Specifically, here we will be dealing with the relation between the Mellin moments of 
the coefficients appearing in the expansion 

A(N) = f 1 dxx N ~ 1 A(x) 
Jo 

B(N) = I dxx N - 1 B(x) 
Jo 
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(1.73) 

and those of the distributions 

A T q( ± \N,Q 2 )= f 1 dxx N - 1 A T q( ± \x,Q 2 )). (1.74) 
Jo 

For this purpose we recall that the general (nonsinglet) solution to NLO for the latter 
moments is given by 



-2A T P$(N)/p 



A T q ± (N,Q 2 ) = K (QiQ\N)^-^j A T q ± (N,Q 2 ) 
with the input distributions A T q±(Ql) at the input scale Q . We also have set 



K(Ql Q\ N) = l + ^Ql)- - S m W tP V {n) _ £La t P^(N) 

7rpo ^Po 



(1.75) 



In the expressions above we have introduced the corresponding moments for the LO and 
NLO kernels (A T P$' N , A T P$£). 

The relation between the moments of the coefficients of the nonsinglet x-space expan- 
sion and those of the parton distributions at any Q, as expressed by eq. (1.75) can be 
easily written down 



A n (N) + a s B n {N) = A T q ± (N, Q 2 )K(Q , Q, N) i^lM^j . (1.76) 

As a check of this expression, notice that the initial condition is easily obtained from 
(1.76) setting Q — > Qo,n — > 0, thereby obtaining 

A» S (N) + a s B» s (N) = A T q ± (N, Q 2 ), (1.77) 

which can be solved with A$ S (N) = A T q ± (N,Q 2 ) and Bg s (N) = 0. 

It is then evident that the expansion (1.60) involves a resummation of the logarithmic 
contributions, as shown in eq. (1.76). 

In the singlet sector we can work out a similar relation both to LO 

^ Ait) Ait)" 

with 

* = AiS pm(N) -^ 
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e 2 = -J— (-pW^ + Axl 

A2 — A\ v 



Ai, 2 = \ (P£\N) + PW(JV) ± sj (P£\N) - P!S\N)) 2 + 4P£\N)p£\N)\l.79) 



and to NLO 



where 



A n (N) + a s B n (N) = xi Mr 1 )" + X2 ("^V , (1-80) 



a (— 2 e 2 Rex 
Xi = e i + tt -5-eiRei + 



2n\(3 A 1 -A 2 -/3 /2 / 

a /-2 exRea 
X2 = e 2 + — — e 2 Re 2 + 



with 



2vr V/?o A 2 -A 1 -/3 /2 / 



R = pW(jV) - Ap(o)(jv). (1.82) 

We remark that A n (N) and B n (N), P (0) (iV), pW(JV), in this case, are all 2-by-2 singlet 
matrices. 

1.8 Initial conditions 

As input distributions in the unpolarized case, we have used the models of Ref.[71], valid 
to NLO in the MS scheme at a scale Ql = 0.40 GeV 2 

x(u-u)(x,Ql) = 0.632a; a43 (l -x) 3m (l + 18.2a;) 

x(d-d)(x,Q 2 ) = 0.624:(l-x) 10 x{u-u)(x,Ql) 

x(d-u)(x,Ql) = 0.20a; - 43 (l - x) 12A (l - 13.3^ + 60.0a;) 

x(u + d)(x,Ql) = 1.24a; - 20 (1 - x) 8 - 5 (1-2.3^ + 5. 7x) 

xg(x,Q 2 ) = 20.80a; L6 (l - a;) 4 ' 1 (1.83) 

and xqi(x, Ql) = xq^(x, Ql) = for q { = s, c, b, t. 

Following [70], we have related the unpolarized input distribution to the longitudinally 
polarized ones 

xAu(x,Q 2 ) = 1.019x a52 (l -a;) ai2 XM(a;,g 2 ) 
a;Ad(a;,Q 2 ) = -0.669x 0A3 xd(x,Q 2 ) 



(1.81) 
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xAu(x,Ql) = -0.272x - 38 xu(x,Ql) 
xAd(x, Ql) = xAu(x, Ql) 

xAg(x,Ql) = lA19x 1A3 (l-x)°- 15 xg(x,Ql) (1.84) 

and xAqi(x,Ql) = xAqj(x,Ql) = for q; L = s,c,b,t. Being the transversity distribution 
experimentally unknown, following [96], we assume the saturation of Soffer's inequality 

xA T q t (x, Ql) = m(x,Ql)+xA qi (x,Ql) ^ (Lg5) 

These input distributions will be used in the next section in the analysis of the trans- 
verse asymmetries for polarized pp collisions in the Drell-Yan process. 



1.9 Transversity in the Drell-Yan process 

The Drell-Yan mechanism for lepton pair production is, at parton level, described by 
the annihilation of a quark-antiquark pair (qq) into an s-channel (virtual) photon which 
decays into a lepton pair. The process is characterized by a distinct signature, since the 
two leptons of the final state can be more easily identified. As we have discussed in 

the previous sections of this chapter, the study of the transverse spin distributions of the 
proton is an ongoing process which requires more experimental data in the future in order 
to provide us with a clearer understanding of these functions. The natural question to ask 
is: what are the processes mediated at parton level by these particular matrix elements, 
which therefore carry information, over to the final state, on the distribution of transverse 
spin in the proton. In this section we focalize our attention on the phenomenology of the 
transversely polarized distributions h1(x,Q 2 ) (or A T (x, Q 2 ) in other notations) which is 
the missing part in the QCD description of the spin structure of the nucleon at leading 
twist. 

Since h1(x, Q 2 ) is chirally-odd, as discussed by Jaffe, one of the possible way to access 
the measurement of this observable is the Drell-Yan process. 

In fact, in deep inelastic scattering processes h1(x,Q 2 ) is severely suppressed in the 
operator product expansion since it needs a mass insertion in the unitarity diagram to 
appear. So the Drell-Yan process remains up to now the theoretically cleanest way to 
observe Q 2 ). 

In a polarized proton-antiproton collision one can construct an asymmetry Att which 
is proportional to the product of transversity of the proton's quark and the transversity 
of the antiproton's anti-quark, which are equal due to the charge conjugation. The most 
recent experimental proposal of antiproton-proton scattering with polarization has been 
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presented by the PAX collaboration. In the PAX experiment polarized antiprotons are 
produced by spin filtering an internal polarized gas target and scattered off protons at 
intermediate energy. The computation of the related transverse spin asymmetries in Drell- 
Yan, which we are going to discuss below, is performed using the factorization formula 
for the cross section dda = (da^-da n )/2. This is given as a double convolution of 
transversity distributions with the corresponding transversely polarized partonic cross sec- 
tion. In the case of antiproton-proton collision with a di-muon production the expression 
is 



d§U rl /"l _ 

= J2^ 2 q dx! dx 2 A T q(xi,(j, 2 F )A T q(x2,(J? F ) + A T q(xi, ii F )A T q(x 2 , ii f ) 



dMdydcj) 

dA T a 

x — 4 ,, 1.86 

dMdydcp y ' 

where \xp is the factorization scale of the process, and the charge e is quite general, 
since it may encompass both electromagnetic and electroweak effects. S is the center 
of mass initial energy and M represents the invariant mass of the virtual photon. The 
y variable is called rapidity and it is connected with the variables x\, x^ by a scaling 
parameter r = M 2 /S as x\ = \fre y , x° = \fre~ y . is the azimuthal angle of one muon. 

The next-to-leading order a s expression [97] for the hard scattering term is written in 
the modified minimal subtraction MS scheme as follows 



(1),M5 



dAxcr 

dMdydcj) 



x 



4r(xix 2 + r) 



9SM F 2tt x^ix, + x ^ + x° 2 ) { <P) 
\5(x 1 -x 1 )5(x 2 -x 2 ) 



l ln2 (l-*g)(l-s°) vr 2 
4 t 4 



+6(X! - X°) 



(x 2 - x 2 y )+ t(x 2 + xl) 



\n(x2 — X2) 

v X 2 -X°2 



+ 



1 X°2 

— ln — 

x 2 - x 2 x 2 



2[(x 1 -x° 1 )(x 2 -x° 2 )}_ 
[1 - 2] . 



+ 



(xix 2 + x 2 x\) 2 (x\ — Xi)(x 2 — x 2 ) 



} + 



(1.87) 



In this expression the renormalization scale /jLr has been set to coincide with the factor- 
ization scale \i F = M. We recall that a different choice for these two scales would be 
responsible for the generation of additional logs (\og(/i F / both in the evolution and 
in the hard scatterings. This issue will be discussed in more detail in the computation of 
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the Higgs total cross section, where it has more relevance. 

The expression of the LO cross section is quite simple and it reads 

^ = ^-(20)^-,?)^ -A)- (1-88) 

To lowest order (LO) x\ and x® coincide with the momentum fractions carried by the 
incident partons. 

The asymmetry depending on y variable can be constructed 



' Jl) ' " dM fi* d<f> da/dMdyd<j) ' 



A measurement of the asymmetry with a sufficient accuracy can be done in the PAX 
experiment in the dilepton mass region below the J/\l/ threshold and a systematic study 
at leading order has be done in this last two years [51, 9]. We have numerically analyzed 
[15] a region very close to the J/\l/ resonance Mw4 GeV 2 to NLO since could be crucial 
in order to enhance the cross section. In fact asymmetries would be very difficult to 
measure in a region characterized by a fast falling cross section. Near a resonance the 
possibility of a successfull measurement is sharply enhanced. The cross section for dilepton 
production increases by almost 2 orders of magnitude in going from M = 4 to M = 3 GeV, 
since this cross section involves unknown quantities related to qq-J/^ coupling. However, 
independently of these unknown quantities, the qq-J/^ coupling is a vector one, with the 
same spinor and Lorentz structure as qq-j* coupling. These unknown quantities cancel in 
the ratio giving Att, while the helicity structure remains, so that the asymmetry is given 
by 

^(x 1; M 2 K(a; 2 ,M 2 ) 

A TT ~ a TT — ( tj^— ( tj^- , (1.90) 

u(xi, M 2 )u(x 2 , M l ) 

where cltt is the double spin asymmetry of the elementary QED process qq — > / + /~ 

sin 2 6 
1 + cos 2 

This substantially enhances the sensitivity of the PAX experiment to Att and the amount 
of direct information achievable on /i"(xi, M 2 )h^(x2, M 2 ). In order to evolve the transver- 
sity we have used the longitudinal bound on the inputs 

hl(x,Q 2 ) = Aq(x,Q 2 ) (1.92) 
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Figure 1.1: NLO asymmetry calculated at M = 4 GeV with MRST inputs. 



and we performed the evolution at LO and NLO starting by the same initial point Ql = 1 
GeV 2 for the GRV's and Q 2 = 0.4 GeV 2 for the set of distributions termed MRST's. 
We show below some plots of the asymmetry vs the rapidity calculated at NLO using 
different input parton distribution functions and at different values of S. In fig. (1.1) we 
present NLO results for the transverse double spin asymmetries for an invariant mass of 
the lepton pair equal to 4 GeV and in fig. (1.2) we plot the unpolarized cross section at 
the same energy. The asymmetry is about 15 percent, considering a GRV input model 
for the transverse spin distributions, with a sizeable cross section. We also show results 
for the integrated asymmetries in figs. (1.3) and (1.4), which grow up to 35-40 percent, 
and are therefore sizeable. 



1.10 Documentation of the Code 



In this section we describe the variables, the parameters and the functions introduced in 
the numerical implementation of the program. The notation that we use is the standard 
one adopted by Computer Physics Communication for the documentation of the code. 
This section consists of a list of all the functions and variables defined in the program, 
the output files, and some comments concerning the performance of the implementation. 
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Figure 1.2: NLO cross section calculated at M = 4 GeV with MSRT inputs. 

1.10.1 Names of the input parameters, variables and of the out- 
put files 



1.10.1 


1 Notations 







gluons, g 


g 


1-6 


quarks, q i7 sorted by their mass values(w, d, s, c, b, t) 


u,d, s , c ,b,t 


7-12 


antiquarks, ql 


au , ad , as , ac , ab , at 


13-18 




um , dm , sm , cm , bm , tm 


19-24 


Xi (unpolarized and longitudinally polarized cases) 


Cu,Cd,Cs .Cc^b.Ct 




q^ (transversely polarized case) 


Cu,Cd,Cs ,Cc,Cb,Ct 


25 


q (+) 


qp 



1.10.1.2 Input parameters and variables 

process unpolarized 

1 longitudinally polarized 

2 transversely polarized 
spacing 1 linear 

2 logarithmic 

GRID_PTS Number of points in the grid 

MGP Number of Gaussian points, uq 

ITERATIONS Number of terms in the sum (1.60) 

extension Extension of the output files 
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Figure 1.3: 
inputs. 



NLO Asymmetry calculated integrating M between 2 and 3 GeV with GRV 



step 


grid step (linear spacing case) 


lstep 


step in log 10 x (logarithmic spacing case) 


X[i] 


i-th grid point, Xi 


XG[i] [j] 


j-th Gaussian abscissa in the range [x^, 1], 


WG[i] [j] 


j-th Gaussian weights in the range [xi, 1], 


nf, Nf 


number of active flavors, n/ 


n_evol 


progressive number of the evolution step is rif — 3 


Q[i] 


values of Q in the corresponding grid 


lambda [i] 


' wnere i = nf — 3 


A[i] [j] [k] 


coefficient Aj(xk) for the distribution with index % 


B[i] [j] [k] 


coefficient Bj(xk) for the distribution with index i 


betaO 


A) 


betal 


A 


alphal 


&s(Qin), where Qi n is the lower Q of the evolution step 


alpha2 


® s (Qfin), where <5/m is the higher Q of the evolution step 



1.10.1.3 Output files 

The generic name of an output file is: Xni.ext, where 



X is U in the unpolarized case, L in the longitudinally polarized case and T in the trans- 
versely polarized case; 
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Figure 1.4: NLO Asymmetry calculated integrating M between 4 and 7 GeV with GRV 
inputs. 

n is a progressive number that indicates the scale Q 2 at a given stage: n = refers to 
the initial scale, the highest value of n refers to the final scale and the intermediate 
values refer to the quarks production thresholds (1 for charm, 2 for bottom and 3 
for top); 

i is the identifier of the distribution, reported in the third column of the table in subsub- 
section 1.10.1.1; 

ext is an extension chosen by the user. 

1.10.2 Description of the program 
1.10.2.1 Main program 

At run time, the program asks the user to select a linear or a logarithmic spacing for the 
x-axis. The logarithmic spacing is useful in order to analyze the small-x behavior. Then 
the program stores as external variables the grid points Xi and, for each of them, calls the 
function gauleg which computes the Gaussian points Xy and weights Wij corresponding 
to the integration range [xi, 1], with < j < no — 1- After that, the user is asked to enter 
the type of process, the final value of Q and an extension for the names of the output 
files. At this point the program computes the initial values of the parton distributions for 
gluons, up, down, anti-up and anti-down (see Section 1.8) at the grid points and stores 
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them in the arrays A[i] [0] [k] (see (1.64)), setting to zero the initial distributions of the 
heavier quarks. 

The evolution is done in the various regions of the evolutions, all characterized by 
a specific flavour number. Each new flavour comes into play only when the scale Q 
reaches the corresponding quark mass. In that case n/ is increased by 1 everywhere in 
the program. The recurrence relations (1.62) and (1.63) are then solved iteratively for 
both the nonsinglet and the singlet sector, and at the end of each energy step the evolved 
distributions are reconstructed via the relation (1.60). The distributions computed in this 
way become the initial conditions for the subsequent step. The numerical values of the 
distributions at the end of each energy step are printed to files. 

1.10.2.2 Function writef ile 

void writef ile (double *A,char *f ilename) ; 

This function creates a file, whose name is contained in the string *f ilename, with 
an output characterized by two columns of data: the left column contains all the values 
of the grid points Xi and the right one the corresponding values of the array A [i] . 

1.10.2.3 Function alpha_s 

double alpha_s (double Q, double betaO, double betal, double lambda); 

Given the energy scale Q, the first two terms of the perturbative expansion of the 
/5-function betaO and betal and the value lambda of A^4, alpha_s returns the two-loop 
running of the coupling constant, using the formula (1.5). 

1.10.2.4 Function gauleg 

void gauleg(double xl, double x2, double x[], double w[],int n) ; 

This function is taken from [104] with just some minor changes. Given the lower and 
upper limits of integration xl and x2, and given n, gauleg returns arrays x[0, . . . ,n-l] 
and w [0 , . . . ,n-l] of length n, containing the abscissas and weights of the Gauss-Legendre 
n-point quadrature formula. 

1.10.2.5 Function interp 

double interp (double *A, double x) ; 

Given an array A, representing a function known only at the grid points, and a number 
x in the interval [0, 1], interp returns the linear interpolation of the function at the point 
x. 
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1.10.2.6 Function IntGL 

double IntGL(int i, double kernel(double z), double *A) ; 

Given an integer i (corresponding to a grid point Xi), a one variable function kernel (z) 
and an array A, representing a function g(x) known at the grid points, IntGL returns the 
result of the integral 

I —kernel ( — ] g(xi), (1.93) 
computed by the Gauss-Legendre technique. 



1.10.2.7 Function IntPD 

double IntGL (int i, double *A) ; 

Given an integer i, to which it corresponds a grid point Xi, and an array A, representing 
a function f(x) known at the grid points, IntGL returns the result of the convolution 

1 ^ tf \ f 1 d yyf(y) - x if( x i) . m * n nA ^ 

v- (8) /(Xi) = / + /(a: i )log(l-a: i ), 1.94 

(1 -Xi)+ Jxi y y-Xi 

computed by the Gauss-Legendre technique. 



1.10.2.8 Function S2 

double S2 (double z) ; 

This function evaluates the Spence function ^(z) using the expansion 

S 2 (z) = logzlog(l -z)-- (logz) 2 + - + E (!- 95 ) 

4 iz n =l 71 

arrested at the 50th order. 



1.10.2.9 Function fact 

double fact (int n) ; 

This function returns the factorial n\ 



1.10.2.10 Initial distributions 

double xuv (double x) ; 
double xdv (double x) ; 
double xdbmub (double x) ; 
double xubpdb (double x) ; 
double xg (double x) ; 
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double xu (double x) ; 
double xubar (double x) ; 
double xd (double x) ; 
double xdbar (double x) ; 
double xDg (double x) ; 
double xDu (double x) ; 
double xDubar (double x) ; 
double xDd (double x) ; 
double xDdbar (double x) ; 

Given the Bjorken variable x, these functions return the initial distributions at the 
input scale (see Section 1.8). 

1.10.2.11 Regular part of the kernels 

double PONS (double z) ; 
double P0qq(double z) ; 
double POqg (double z) ; 
double P0gq(double z) ; 
double POgg (double z) ; 
double PlNSm(double z) ; 
double PINSp (double z) ; 
double Plqq(double z) ; 
double Plqg (double z) ; 
double Plgq(double z) ; 
double Plgg (double z) ; 
double DPONS (double z) ; 
double DPOqq (double z) ; 
double DPOqg (double z) ; 
double DPOgq (double z) ; 
double DPOgg (double z) ; 
double DP lMSm (double z) ; 
double DPlMSp(double z) ; 
double DPlqq(double z) ; 
double DPlqg(double z) ; 
double DPlgq(double z) ; 
double DPlgg (double z) ; 
double tPO (double z) ; 
double tPlm(double z) ; 
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double tPlp (double z) ; 

Given the Bjorken variable z, these functions return the part of the Altarelli-Parisi 
kernels that does not contain singularities. 

1.10.3 Running the code 

In the plots shown in this chapter we have divided the interval [0, 1] of the Bjorken variable 
x in 500 subintervals (GRID_PTS=501), 30 Gaussian points (NGP=1), and we have retained 
10 terms in the sum (1.60) (ITERATI0NS=10). In the figures 1.11 and 1.16 the flag spacing 
has been set to 2, in order to have a logarithmically spaced grid. This feature turns useful 
if one intends to analyze the small-x behavior. We have tested our implementation in a 
detailed study of Soffer's inequality up to NLO [29]. 
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Figure 1.5: Evolution of the unpolarized quark up distribution xu versus x at various Q 
values. 
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Figure 1.6: Evolution of xd versus x at various Q values. 
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Figure 1.7: Evolution of xs = xs versus x at various Q values. 
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Figure 1.8: Evolution of the unpolarized antiquark up distribution xu versus x at various 
Q values. 

1.11 Unpolarized kernels 

The reference for the unpolarized kernels is [63], rearranged for our purposes. We remind 
that the plus distribution is defined by 



o (1 — x)+ Jo 



1 — X 



(1.96) 



and the Spence function is 



S 2 (x) = -2Li 2 (-ar) - 2 logrr log(l + a;) + -log 2 :r- —, 



(1.97) 



where the dilogarithm is defined by 



Li2W= /°Ml^) dt . 

J x t 



(1.98) 
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Figure 1.9: Evolution of the unpolarized gluon distribution xg versus x at various Q 
values. A 




P$l-(x) = P$l + (z) = P£\x)=C F 



P£\x) = 2T f [x 2 + (l-xy 



X 



P^{x)=2N c 



+--2 + x(l-x) 



[1 — x) + X 



+ ^6(1 -x) 



(1.99) 

(1.100) 
(1.101) 

(1.102) 



Pns- ( x ) 



Cf 
18 



162C F (x - 1) + 4T f (llx - 1) + N C (S9 - 223a; + 3vr 2 (l + a:)) 
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X 



Figure 1.11: As in figure (1.10), but now the x-axis is in logarithmic scale, to show the 
small- a; behavior. 
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2{x-l) 



log 2 a; 
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C F (2C F -N c )(l + x 2 ) 



1 + x 



loga;log(l — x) 

S 2 (x) 



9l 

9 

C F 
{ ~2 



N c (3n 2 - 67) + 207) 
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log a; 



(1.103) 
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Figure 1.13: Evolution of xAs = xAs versus x at various Q values. 
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Figure 1.14: Evolution of the longitudinally polarized antiquark up distribution xAu 
versus x at various Q values. 



Pns+ ( x ) 



18 



18C F (x - 1) + 4T f (llx - 1) + N c (17 - 151x + 3tt 2 (1 + x)) 



+ 
+ 
+ 
+ 



C F [6CfQ. + 2x) - (11N C - 47)) (1 + x 2 )} 

6(x-l) 
C F [C F - N c - (C F + N c )x 2 } 



log a; 



2cm + x 2 



2(x - 1) 

logxlog(l — x) 



log 2 X 



x — 1 

C F (2C F -N c )(l + x 2 ) 



1+x 



S 2 (x) 



C f 



+ 



9 

C f 



N c (3n 2 - 67) + 20T f 



l-x) + 

[N c (51 + 44tt 2 - 216C(3)) - 47) (3 + 4tt 2 ) 



46 



1.11. Unpolarized kernels 



0.25 



0.15 




0.05 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 1.15: Evolution of the longitudinally polarized gluon distribution xAg versus x at 
various Q values. 
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Figure 1.16: As in figure (1.15), but now the x-axis is in logarithmic scale. 
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(1.104) 
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Figure 1.17: Evolution of the transversely polarized quark up distribution xA F u versus x 
at various Q values. 
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Figure 1.18: Evolution of xAtcI versus x at various Q values. 
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(1.105) 



pW(x) = |^[T / (3C F a;(42 - 87a; + 60a; 2 -7r 2 (2 + 4(a;-l)a;)) 

+iV c (40 + a;(450:r - 36 - 436a; 2 + tt 2 (3 + 6(x - l)x))))] } 

+ |^- [6N C + 8N c x(6 + llx) + 3C F (3 - Ax + 8x 2 )] J logx 

+ {8(C F - N c )T f (l - x)x} log(l - x) 

+ [T } [C f (1 -2x + Ax 2 ) - N c (3 + 2x(3 + x))] } log 2 x 

+ {2(C F - N c )T f [1 + 2(x - l)x}} log 2 (l - x) 

- {AC F T f [l + 2(x-l)x]} logxlog(l - x) 

+ {2N c T f [1 + 2x(l + x)}} S 2 (x) (1.106) 
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Figure 1.19: Evolution of the transversely polarized antiquark up distribution xA F u versus 
x at various Q values. 
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X 



S 2 (x) 



(1.107) 



Pff(x) = {^[24C F T f (x -l)(x(ll + 5x) -l)+4N c T f (x(29 + x(23x - 19)) - 23) 

+N% (67r 2 (a;(2 + (x - l)x) - 1) - x(25 + 109a;))] } 
+ [ Ag [11(1 - 4:r):r - 25] - 4N c T f (l + x) - 6C F T f (3 + 5s) j ^ 

f 2C F T / a:(a: 2 - 1) + iV 2 [1 + x (2 + g(3 + (x - 6)z))] 1 2 ^ 
I (1 - x)x J ° g X 

[47V 2 [l + (x-l)x] 2 l , , n . 
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2iV 2 (1 + x + x 2 ) 



2^2' 



x(l + x) 



S 2 (x) 



+ 



N c 

IT 

N c 
IT 



N c {3n 2 - 67) + 207) 

J J (l-a:)+ 
[7V C (8 + 9C(3)) - 47}] - C F 7)} 5(1 - x) 



(1.108) 



1.12 Longitudinally polarized kernels 



The reference for the longitudinally polarized kernels is [114]. 



,(o) 



2 3_, 
(1^)7 ~ l ~ x+ 2 ^ ) 



AP£\x)=2N c 



AP£\x) = 2T,(2x-l) 

ApV>(z) = C F (2-z) 
1 



(1-x) 



-2x + l 



+ ^S(l-x) 
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AP«_(*)=P« + (:r) 
AP£ + (x)=?S_(x) 

162C F (a; - 1) + 8T)(4 + z) + N c (89 - 223a; + 3tt 2 (1 + x 

N c (x 2 - 23) - QC F (4x 2 -2x-5) 
+87) (2 + x{5x - 6))]}logx 



(1.109) 

(1.110) 
(1.111) 

(1.112) 

(1.113) 
(1.114) 



C F 



6(x- 1) 



+ 



' CV [C F - N c + 47) - (Cf + N c + 47))a; 2 
2(x-l) 



• log 2 x 



+ {^^!)}log,log(l-,) 

-{ CF(2CF l+f 1+x2) }^ 

-{^[iVc(37r 2 - 67) + 20T / ]}^ 
I 72 



;r) 



+ 



+ 



N c (51 + 44tt 2 - 216C(3)) - 47) (3 + 4tt 2 ) 
+9C F (3 - 4tt 2 + 48C(3))] } 5(1 - x) 



(1.115) 
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C F (vr 2 (2 - Ax) - 66 + Six) + N c (72 - 66a; + vr 2 (2x - 1) 



+ {Tf [2N C (1 + 8x) - 9C F ]}logx 

+ {8(N C - C F )T f (x - 1)} log(l - x) 

+ {Tf [C F (2x - 1) - 3iV c (l + 2x)]} log 2 x 

+ {2(C F - N c )Tf(2x - 1)} log 2 (l - x) 

+ {AC F T f (l -2x)} log a; log(l - x) 

+ {2N c T f (l + 2x)}S 2 (x) 



(1.116) 



+ 



+ 



{if [ 9C ^(8x - 17) - 87>(4 + x) + N c (82 + 3tt 2 (x - 2) + 70x)] J 
[iV c (8 - 26x) + C F (x - 4)]} logo: 
|^ [47)(x - 2) - 3C F (2 + x) + iV c (10 + x)]} log(l - x) 

+ [3iV c (2 + x) - C F (x - 2)]} log 2 x 

+ {CV(CV - AT c )(x - 2)} log 2 (l - x) 
+ {2C F N c (x - 2)} logxlog(l - x) 

~{C F N c (2 + x)}S 2 (x) (1.117) 



180C F T/(x - 1) + 8N c T f (19x - 4) + N c (6tt 2 (1 + 2x) - 305 - 97x)] | 

+ |1 [AT 2 (29 - 67x) + 6CV7)(x - 5) - 4JV C 7)(1 + x)] } logx 
' N c (2x 2 + x - 4) - 2C F T)(x 2 - 1) " 



+ 



x- 1 



log X 



'4iV 2 x(2x- 1)1 , , „ 
f <! — \ logxlog(l - x) 



x- 1 
'2iV 2 x(l + 2x) 



S 2 (x) 



N c 



1 + x 

"iVc(37r 2 - 67) + 20T 



+ 



(l-x) + 

1 1 \nI ( 8 + 9C(3)) - 3C F T f - AN c Tf] } 5(1 - x) 



(1.118) 



1.13 Transversely polarized kernels 

The reference for the transversely polarized kernels is [115]. 
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A T p£l(x) = A T P%> + (x) = C F 



(0) 
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1 + x 

N c (3n 2 - 67) + 207)' , , , 

AT c (5i + 44^ _ 216C(3)) - 47)(3 + 4tt 2 ) 
+9C F (3 - 4tt 2 + 48C(3))] } 5(1 - x) 



(1.120) 
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3(x- 1) 

' 4C Mi i n ^ 
> log a; log (1 — 

2C F (iV c -2C F )a;' 



1 + x 



+ 



9 

Cf 
72 



iV c; (37r 2 - 67) + 20T / 

J JJ (l-x)+ 

7V C (51 + 44vr 2 - 216C(3)) - 47) (3 + 4vr 2 ) 
+9C F (3 - 4tt 2 + 48C(3))] } 5(1 - x) 



(1.121) 



1.14 Conclusions 



We have illustrated and documented a method for solving in a rather fast way the NLO 
evolution equations for the parton distributions at leading twist. The advantages of the 
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method compared to other implementations based on the inversion of the Mellin moments, 
as usually done in the case of QCD, are rather evident. We have also shown how Rossi's 
ansatz, originally formulated in the case of the photon structure function, relates to the 
solution of DGLAP equations formulated in terms of moments. The running time of the 
implementation is truly modest, even for a large number of iterations, and allows to get 
a very good accuracy. 



Chapter 2 



NNLO extension of the solution of the 
Renormalization Group Equations 

2.1 Introduction to the Chapter 

We have seen in Chapter 1 that renormalization group equations play a key role in the 
computation of cross sections in hadronic collisions. We have also seen that the depen- 
dence on the factorization scale of the factorization formula is largely reduced when higher 
order corrections in 0:5 are kept into account. As we have already discussed, the use of 
efficient algorithms for the solution of these equations is of remarkable relevance in or- 
der to efficiently compare the theory with the experiments. As we increase the collision 
energy in a scattering process, the factorization scale has a wider interval allowed for its 
variation. It is of some help to summarize this issue in a simple but rather clarifying way. 

Consider a proton proton collision at LHC energy (these can cover a range between 
few TeV's and 14 TeV) and let us assume that we use the parton model description of the 
process based on the factorization formula. As we have mentioned, Q, the factorization 
scale (sometimes also called "///") is a fraction of the total energy available in the center 
of mass frame of the two colliding beams. In order to decide on the most accurate 
value for \if so to come up with a prediction, let's say, for a multi jet cross section, 
we have to look at the average pr of the final state jet and use that energy value as 
an indication for the underlying value of Of course this procedure is approximate 
and requires a study of the behaviour of the final result for the cross section in terms of 
various choices of ///. A drastic reduction on the dependence of the result on the value 
of /I f is obtained if we include higher order corrections. As we have already stressed 
before, both the hard scatterings and the evolution of the parton distributions have to 
be accurate enough for this reduced sensitivity to emerge. In general the computation of 
hard scatterings, which are process dependent, are laborious enough to be limited to a 
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specific set of "golden plated modes", such as Drell-Yan lepton pair production and few 
more, and require difficult technical studies by various research groups to be performed. 
For instance, in the case of the Higgs total cross section and of the cross section for the 
production of a Higgs particle with an associated gauge boson at the LHC, considerable 
progress has been done in the last few years. On the other hand the computation of the 
kernels of the renormalization group equations for the parton distributions has also to be 
known at the same level of accuracy as the hard scatterings. Therefore it is imperative, 
especially for the detection of the Higgs boson at LHC energy, to be able to move from 
the NLO to the NNLO case in the study of some "golden plated modes" involving this 
particle. It is unlikely that in the near future most of the NLO hard scatterings will be 
extended to include their NNLO corrections, but, as we have mentioned, in the case of 
the Higgs this process is under intensive investigation. The computation of the NNLO 
evolution kernels, which are process independent, has been completed quite recently and 
has set a landmark on the applicability of perturbative QCD to hadron colliders. The 
numerical study of these kernels is a nontrivial task and the accurate solution of the 
associated RGE's is also nontrivial. 

In this second chapter our aim is to extend up to NNLO the results of chapter 1 
using a similar methodology in order to solve these equations. At this time our numerical 
implementation is the second code which is able to evolve parton distributions up to this 
level of accuracy. It is obvious that most of the discussion in this chapter overlaps with 
the methods of chapter 1, but with some key differences that we will underline. The new 
implementation contains new features which are not present in our previous discussion 
since the new NNLO kernels are far more complex than those known at NLO. These issues 
are discussed in some detail in this and in the next chapter. 

After the submission of this thesis, the preliminar work presented in this chapter has 
been completed and published in [30]. 



2.2 Definitions and Conventions 

In this section we present our definitions and conventions, which are similar to those of 
chapter 1 with some important modifications. We introduce the 31oop evolution of the 
coupling [36] 



a s (Q 2 ) 



An 



i _AlogL + 1 



P 2 o L ffil? 



^(log 2 L-logL-l) +/ 3 2 + C>( 



+ OI.Z3 



(2.1) 



where 



(2.2) 
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bb 



and the beta function is defined by 



da s {Q 2 ) 

p(as) = dW (2 - 3) 

and its three-loop expansion [113] is 
where 

Pi = y n c - Y Ncnf - 2Cpn ^ ^ 

= + 2C F T f - 2 -fc F N c T f - ±§^ T/ + fc F Tf + l -^N c T] (2.7) 

and 

NZ-1 4 1 

^c = 3, C F = ^— = -, = T R rif = -rtf, (2.8) 

where Nc is the number of colors, n/ is the number of active flavors, that is fixed by 
the number of quarks with m q < Q; A-j^ is calculated using the known value of a s (mz) 
and imposing the continuity of a s at the quark masses thresholds. We recall that the 
perturbative expansion of the kernels now includes the NNLo contributions is 

P(x,a s ) = (g) pW(x) + V>(s) + (^) 3p(2) (-) + • • • • (2-9) 
whose specific form can be found in the original literature [100, 117]. 



We solve Eq. (1.13) directly in x-space, assuming a solution of the form 

+ {a , m yf;£M log ^jm (2 , 0) 



n=0 



Us{Ql) 



for each parton distribution /, where Qo defines the initial evolution scale. 

As in chapter 1, also in this case we derive the following recursion relations for the 
coefficients A n , B n and C n (see the frame below for details) 

A n+1 (x) = -^-P (0 \x) ® A n (x), (2.11) 
Po 
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B n+l {x) = -B n {x) - -^-A n+1 (x) - ^P (0 \x) ® B n (x) - -^-P (1 \x) <g> A n (x), (2.12) 



C n+1 (x) = -2C n (x) - ^B n (x) - ^B n+1 {x) - J^-A n+1 (x) 

-^P (0) (x) ® C n (x) - -^-P {1 \x) ® B n (x) 
Po ttPo 

l —P i - 2 \x)®A n {x). (2.13) 



Derivation of the recursion relations and renormalization scale depen- 
dence. 

We introduce the shortcut notation 

L ^>^ a ^m (2 - 14) 

and, making use of the beta function definition (2.3), we compute its derivative 
dL(Q 2 ) a s {Ql) d a s (Q 2 ) 1 da s (Q 2 ) f3{a s ) 



dlogQ 2 a s (Q 2 )d\ogQ 2 a s (Q 2 ) a s (Q 2 ) dlogQ 2 a s (Q 2 ) 



(2.15) 



Inserting our ansatz (2.10) for the solution into the DGLAP equation (1.13) we get 
for the LHS 

y l Mx) nL n-lPM + a B n (x) nL n-lPM 

~y \ n\ a s n\ a s 



+a, ; — nL 

n\ a s J 

+ E IPM^L" + 2a.P( a .)^L») . (2.16) 

Note that the first sum starts at n — 1, because the n = term in (2.10) does not 
have Q 2 dependence. Tranforming n — > n — 1 in the first sum, using the three-loop 
expansion of the beta function (2.4) and neglecting all terms of order a 4 s or more, the 
previous formula becomes 



a, — — — to, 



n=0 



4tt 16tt 2 s 64tt 3 



, ^n+l(g) r n / A) 2 Pi „A , Cn+l(g) rw / 
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, B n{x) (_Po 2 

n\ \ An s 16n 2 ~ s 



A, 

An 



(2.17) 



Using the kernel expansion (2.9), we get for the RHS 



+ 



a: 



S7T 



An 2 



P (2) ® A r 



P (1) <g) B r 



Air 2 

«| 
2tt 
of 
2tt 



P (1) <g> Ar, 



(x) 



X) 



[p (0) <g> B 
(x) + f- [PW ® C n ] (x) 



(2.18) 



Equating (2.17) and (2.18) term by term and grouping the terms proportional respec- 
tively to a s , a 2 and a z s we get the three desired recursion relations (2.11), (2.12) and 
(2.13). 



Setting Q = Qq in (2.10) we get 

f(x,Q 2 ) = A (x)+a s (Q 2 )B (x) + (a s (Q 2 )f C (x). (2.19) 

Any boundary condition satisfying (2.19) can be chosen at the lowest scale Qo and in our 
case we choose 

B (x) = C (x) = 0, f{x, Ql) = A (x). (2.20) 

What the program does is starting with a parameterized form of the parton distribution 
functions at a low energy scale Q (typically of the order of 1 GeV), computed by some 
specialized groups by fitting experimental data, imposing the boundary condition (2.19) 
for the coefficients A (x), B (x) and Co (a;), computing iteratively A n (x), B n (x) and C n (x) 
up to a certain value of n by the recursion relations (2.11 - 2.13) and then computing the 
sum (2.10). It should be stressed that the solution of the RGE found by this method is 
characterized by two scales, the final evolution scale Q. However, the hard scatterings can 
be renormalized at a scale different from \ip and this clearly introduces a new scale in the 
hadronic cross section which also is an artifact of the perturbative expansion. Keeping 
track of this new scale dependence requires a careful analysis of the perturbative expansion 
of the kernels, as we are going to illustrate below. 

In order to study the renormalization scale dependence of the NNLO kernel it is 
convenient to solve the RGE for the coupling constant thereby expressing the coupling at 
a scale \ip in terms of the coupling at a different point, /Ir, which is the renormalization 
point of the theory, which is also arbitrary 
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1 _ 1 + ^ (ft 



b\ In 



t 1 + M s (/4)] 



(2.21) 



where a s (/i 2 ) = a s (/i 2 ) / (47r). It is possible to introduce a renormalization scale dependence 
in a s through a Taylor expansion of a s (n F ) in terms of a. s (n R ) 



a s (/4) = a 3 (fj? R ) - 



47T (4^)^ 



(2.22) 



where the dependence is shifted into the factor L = ]n(fj, F / fj, R ) Then, a more handy ex- 
pression for the a s coupling obtained by an inverse power of L\ = \n(fi 2 R / Aq CD ) expansion 
can be used for a s (/i R ). Up to NNLO this gives 



47T 



(PoLa) 2 



6i In La + \ [fe 2 (in 2 L A - lnL A - l) + b 2 



(2.23) 



and 



Pij(x, /j, 2 f , /j, 2 r ) 



, - : _ a s (fi R ) p (0) 
4tt 



(4tt) 2 



(4tt) 3 
I) 



P^(x) - 2p LP$\x) - (fhL - P 2 L 2 ) P$\x) 
P,f (x) - SfoLP§\x) - h(3iL - 3(3 2 L 2 ) P«(x 



(4vr) 4 

- (&L - 5/2 ft A>L 2 + P 3 L 3 ) Pjp(x) 



(2.24) 



These new linear combinations of kernels carry an explicit dependence on /i F and are 
suitable for the study of the dependence of the hadronic cross section on /i R . As we 
have mentioned, for a complete picture to emerge, the NNLO corrections to the hard 
scatterings have to be computed in generality, with ji R and /xp held distinct. 
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2.3 Nonsinglet and singlet structure of the kernels up 
to NNLO 

Let us first introduce the notations 

n f 

5?° = ft ±9,, g (±) =E% (±) - ( 2 - 2 5) 

i=i 

The general structure of the nonsinglet splitting functions is given by 

Pm k = P~m k = faPZ + P q l (2-26) 

P q l q k = PqiQk = SikPqq + Pqq- (2-27) 

This leads to three independently evolving types of nonsinglet distributions: the evolution 
of the flavor asymmetries 

q%k = # } - f?° (2-28) 
and of linear combinations thereof is governed by 

P± s = pv ± P%. (2.29) 

The sum of the valence distributions of all flavors evolves with 

PNS = Pqq - Pqq + % (^qq ~ Pqq) = % + P N S - (2-30) 

The quark-quark splitting function P qq can be expressed as 

P qq = P£s + n f (P q s q + PQ = P£ s + P ps . (2.31) 

The nonsinglet contribution dominates Eq. (2.31) at large x, where the pure singlet term 
P ps = Pqq + Pgq is very small. At small x, on the other hand, the latter contribution takes 
over, as xP ps does not vanish for x — >• 0, unlike xP^ s . The gluon-quark and quark-gluon 
entries in Eq. (1.18) are given by 

P q9 = n f P m , (2.32) 

Pgq = P gqi (2.33) 

in terms of the flavor-independent splitting functions P qi9 = P qi9 and P gqi = P g q i . With 
the exception of the first order part of P qg , neither of the quantities xP qg , xP gq and xP gg 
vanishes for x — > 0. 

In the expansion in powers of a s (2.9), the flavor-diagonal (valence) quantity Pq q is 
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of order a s , while P^ and the flavor-independent (sea) contributions P qq and P^ are of 
order a 2 . A non- vanishing difference P£ — P^ occurs for the first time at the third order. 

b o qq qq 

Our next step is to choose a proper basis of nonsinglet distributions that allows us to 
reconstruct, through linear combinations, the distribution of each parton, i.e. the gluon 
distribution g, the quark distributions qi and the antiquark distributions q i7 namely 2n/ + l 
relevant distributions. The singlet evolution gives us 2 distributions, g and q( + \ so we 
need to evolve 2ny — 1 independent nonsinglet distributions. We choose 

1. q(~\ evolving with P^ s ; 

2- <?AT5,ii — — (f° r eac h i such that 2 < i < rif), evolving with P^s'i 
3. q^su = ~ (f° r eac h * sucn that 2 < % < rif), evolving with P^ s . 
We can easily prove that 

1 / n f 

Q ( nL ) • (2-34) 




Choosing i — 1 in (2.34), we compute q^ from the evolved nonsinglets of type 1 and 2 
and q[ + ^ from the evolved singlet q^ and nonsinglet of type 3. Then from the nonsinglets 
2 and 3 we compute respectively q\ ^ and q\ + ^ for each % such that 2 < i < n/, and finally 
% and <?j. 

Life can be made easier going down from NNLO to NLO, as we have P q s q ^ = Pq^ ■ 
This implies (see Eq. (2.30)) that P^'^ = Pn's^i i- e - the nonsinglets q^ and qi?s,ik ev °l ve 
with the same kernel, and the same does each linear combination thereof, in particular 
q\ ^ for each flavor i. The basis of 2nj — 1 nonsinglet distributions that we choose to 
evolve at NLO is 

1. q\ ■* (for each % < rif), evolving with -P/vi^i 

2. q^s,u = — (f° r eacri * suc h ^at 2 < i < rif), evolving with Pjvs^i 

and the same we do at LO, where we have in addition Pps^ = Pns°\ being Pqq^ = 0. 



Remark: the NLO versus the NNLO decomposition 

Prior to moving toward the implementation of this algorithm, which will be discussed 
below, we pause for a moment in order to remark on the differences between the 
singlet /nonsinglet decomposition of the kernels implemented in chapter 1 and those 
discussed in this chapter. 
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In the transverse case there is no coupling between gluons and quarks, so the singlet 
sector (1.18) is missing. This means that A T q^ is a nonsinglet evolving with A T P^ S ; 
but the same kernel evolves also A T q^g ik , so the linear combinations A T q\ + ^ for each 
flavor i. So the 2n/ relevant distributions that we evolve at LO and NLO are 

1. A T q\ -* (for each i < nf), evolving with A T P^ S ; 

2. Arq^ (for each % < rif), evolving with AtP^ s - 



2.4 Harmonic poly logarithms 

We summarize here some key features of the kernel which are relevant for a better under- 
standing of the implementation and refer to the original literature for more details. 

Higher order computations in QCD involve a special type of trascendental functions, 
termed harmonic polylogarithms (HPL). Although their definition is a straightforward 
extension of that of ordinary Spence functions, additional loop integrations in the Feyn- 
man diagram expansion generate multiple integrals of the logarithmic contributions which 
appear at leading order in a s . The study of these functions in physics has been addressed 
in great generality in the last few years and has been quite useful in order to classify 
and study their behaviour. HPL are, in the case of parton distributions and of coefficient 
functions, studied in the standard domain (0, 1), They show branch cuts away from this 
region and can be analytically continued away from it. The NNLO kernel is expressed in 
terms of various polylogarithms. This computation [100, 117] has been performed in the 
A-Mellin space, where the nonsinglet and the singlet anomalous dimensions are related 
to the DGLAP kernel by the well known Mellin transform 

7 (iV)H = - f 1 dxx N - 1 P in \x) . (2.35) 
Jo 

Inverting the above relation, one obtains the expression for the P^ n \x) in the x-space 
that we will use in our implementation. By this procedure the isomorphism between 
harmonic sums of N and HPL is manifest and this can be demonstrated by algebraic pro- 
cedure. Our notation for the harmonic polylogarithms H mii _ jTnw (x) , rrij = 0, ±1 follows 
Ref. [108]. The lowest-weight (k = 1) functions H m {x) are given by 



H (x) = Inx , H±i(x) = =Fln(l=Fx) 



(2.36) 
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The higher-weight (k > 2) functions are recursively defined as 





x 



if m 



1 ,...,m k = 0,...,0 



(2.37) 



mi,...,m k 



[ dzf mi (z)H, 
Jo 




otherwise 



with 



fo(x) 



1 



f±i(x) 



1 



(2.38) 



x 



1 =p x 



#0,...,0,±i,0, 



,0,±i, 



(x) = fl'±(fc+i) )± (i+i),...(a;) • 



(2.39) 



A; 



these identities and many more have been used in the investigation of the kernel. 

2.5 Description of the program 

2.5.1 Main program 

At run time, the program asks the user to choose linear the perturbative order, the 
input model (currently implemented MRST and Alekhin, see Table 2.1), the final value 
of Q and an extension for the output files. Then the program stores as global variables 
the grid points Xi and, for each of them, calls the function gauleg that computes the 
Gaussian abscissas and weights corresponding to the integration range [xj, 1], 
with < j < no — 1. At this point the program stores the initial values of the parton 
distributions at the grid points in the arrays A [i] [0] [k] . 

The evolution is done by energy steps. Each flavor comes into play only when the 
energy in the center of mass frame reaches the corresponding quark mass, so each step 
is determined by the number of flavors involved. The recurrence relations (2.11 - 2.13) 
are then solved iteratively for both the nonsinglet and the singlet sector, and at the end 
of each energy step the evoluted distributions are reconstructed via relation (2.10). The 
distributions computed in this way become the initial conditions for the subsequent step. 
The numerical values of the distributions at the end of each energy step are printed to 
files. 

2.5.2 External files 

Some external files are used by the program. 



1. xpns2e . f and xpij2e . f are Fortran codes by Moch, Vermaseren and Vogt [100, 117] 
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in which are the NNLO kernels are defined. Very few modifications have been done 
to make them compatible with our code. These files need hplog.f to work. 

2. hplog.f is a Fortran code by Gehrmann and Remiddi [66] in which a subroutine 
that computes numerically harmonic polylogarithms up to weight 4 is implemented. 
Harmonic polylogarithms are defined in [108]. 

3. partonww . f is just a merging of the three Fortran codes mrst20011o .f , mrst2001 .f 
and mrstnnlo.f by the MRST group [94, 95] to access to their grids of LO, NLO 
and NNLO parton densities. Very few modifications have been done. 

4. Io2002.dat, alfll9.dat andvnvalfll55.dat are the MRST parton densities grids 
at LO, NLO and NNLO respectively. 

5. a02m.f is the Fortan code by Alekhin [2] to access to his grids of LO, NLO and 
NNLO parton densities. 

6. a02m.pdf s_l_vfn, a02m.pdf s_2_vfn and a02m.pdf s_3_vfn are the Alekhin par- 
ton densities grids in the variable flavor number scheme at LO, NLO and NNLO 
respectively. 

2.5.3 Functions 

2.5.3.1 Function writef ile 

void writef ile (double *A,char *f ilename) ; 

This function creates a file, whose name is contained in the string *f ilename, with 
two columns of data: the left one contains all the grid points Xi and the right one the 
corresponding values A [i] . 

2.5.3.2 Function alpha.s 

double alpha_s(int order, double Q, double lambda); 

(fl f ) 

Given the perturbative order, the energy scale Q and the value lambda of , and 
making use of the values of fa, Pi and /?2, stored as global variables, alpha_s returns the 
running of the coupling constant, using the formula (2.1). 

2.5.3.3 Function gauleg 

void gauleg(double xl,double x2,double x[] , double w[],int n) ; 

This function is taken from [104] with just minor changes. Given the lower and upper 
limits of integration xl and x2, and given n, gauleg returns arrays x[0, . . . ,n-l] and 
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w[0, . . . ,n-l] of lenght n, containing the abscissas and weights of the Gauss-Legendre 
n-point quadrature formula. 

2.5.3.4 Function interp 

double interp (double *A, double x) ; 

Given an array A, representing a function known only at the grid points, and a number 
x in the interval [0, 1], interp returns the linear interpolation of the function at the point 
x. 



2.5.3.5 


Kernels 






double 


PONSCint 


i, double x) 


3 


double 


P0qq(int 


i, double x) 


J 


double 


P0qg(int 


i, double x) 


J 


double 


POgqQnt 


i, double x) 


3 


double 


P0gg(int 


i, double x) 


J 


double 


PlNSm(int 


i, double x 


); 


double 


PlMSp(int 


i, double x 


); 


double 


Plqq(int 


i, double x) 


? 


double 


Plqg(int 


i, double x) 


3 


double 


Plgq(int 


i, double x) 


3 


double 


Plgg(int 


i, double x) 


3 


double 


P2MSm(int 


i, double x 


); 


double 


P2MSp(int 


i, double x 


); 


double 


P2NSv(int 


i, double x 


); 


double 


P2qq(int 


i, double x) 


3 


double 


P2qg(int 


i, double x) 


3 


double 


P2gq(int 


i, double x) 


3 


double 


P2gg(int 


i, double x) 


3 



Given the Bjorken variable x, these functions return, depending on the value of the 
index i 

1. the regular part of the kernel Pi(x); 

2. the plus distribution part of the kernel P2(x); 

3. the delta distribution part of the kernel P%(x), 



as in Eq. (1.22). 
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2.5.3.6 Function convolution 

double convolution(int i, double kernel (int , double) , double *A) ; 

Given an integer i, to which corresponds a grid point Xi, a two variable function 
kernel (i ,x) , representing a kernel, and an array A, representing a function xf(x) = f(x) 
known at the grid points, convolution returns the sum of the three pieces (1.26, 1.30, 
1.27) computed by the Gauss-Legendre technique. 

2.5.3.7 Function RecRel.A 

double RecRel_A (double *A,int k, double P0 (int , double) ) ; 

Given an array A, representing the function A n (x) known at the grid points, an integer 
k (to which corresponds a grid point Xk) and a two variable function P0 (i , x) , representing 
a leading order kernel, RecRel_A returns the RHS of Eq. (2.11) for x = x^. 

2.5.3.8 Function RecRel.B 

double RecRel_B (double *A, double *B,int k, double P0 (int , double ) , 
double PI (int , double) ) ; 

Given two array A and B, representing the functions A n (x) and B n (x) known at the grid 
points, an integer k (to which corresponds a grid point Xk) and two functions P0(i , x) and 
Pl(i,x), representing respectively the LO and NLO part of a kernel, RecRel_B returns 
the RHS of Eq. (2.12) for x = x k . 

2.5.3.9 Function RecRel_C 

double RecRel_C (double *A, double *B, double *C,int k, double P0 (int , double) , 
double PI (int , double) , double P2 (int , double) ) ; 

Given three array A, B and C representing the functions A n (x), B n (x) and C n (x) known 
at the grid points, an integer k (to which corresponds a grid point Xk) and three functions 
P0(i,x), Pl(i,x) and P2(i,x), representing respectively the LO, NLO and NNLO part 
of a kernel, RecRel_C returns the RHS of Eq. (2.13) for x = x k . 

2.5.3.10 Function Li2 

double Li2 (double x) ; 

This function evaluates an approximated value of the dilogarithmic function Li 2 (x) 
using the expansion 

oo n 

Li2{x) = £ - 2 (2.40) 

n=l U 

arrested at the 50th order. 
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2,5.3.11 Function fact 



double fact(int n) ; 

This function returns the factorial n\ 



2.5.4 Distribution indices and identificatives 






gluons, g 


g 


1-6 


quarks, q i7 sorted by increasing mass (u, d, s, c, b, t) 


u,d,s,c,b,t 


7-12 


antiquarks, qj 


au , ad , as , ac , ab , at 


13-18 


<t> 


um , dm , sm , cm , bm , tm 


19-24 




up,dp,sp,cp,bp,tp 


25 


,<"> 


qm 


26-30 


QNS,lii i 7^ 1 


dd, sd , cd,bd ,td 


31 


q (+) 


qp 


32-36 


Q-N~S,lii i 7^ 1 


ds , ss , cs ,bs ,ts 



2.5.5 Input parameters and variables 



GRID_PTS 


Number of points in the grid 


NGP 


Number of Gaussian points, ric 


ITERATIONS 


Number of terms in the sum (2.10) 


extension 


Extension of the output files 


order 


Perturbative order (0=LO, l=NLO, 2=NNLO) 


input 


Input model (1=MRST, 2=Alekhin) 
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X[i] i-th grid point, 

XG[i] [j] j-th Gaussian abscissa in the range [x iy 1], 

WG[i] [j] j-th Gaussian weight in the range [xi, 1], 

nf , Nf number of active flavors, rif 

nf i number of active flavors at the input scale 

Q [i] values of Q between which an evolution step is performed 

lambda [nf] 

MS 

A[i] [j] [k] coefficient Aj(xk) for the distribution with index i 

B[i] [j] [k] coefficient Bj(xk) for the distribution with index i 

C[i] [j] [k] coefficient Cj(xk) for the distribution with index i 

betaO fa 

betal fa 

beta2 p 2 

alphal ai s (Qi n ), where Qi n is the lower Q of the evolution step 

alpha2 &s(Qfm), where Qfi n is the higher Q of the evolution step 

2.5.6 Output files 

The generic name of an output file is: Un_i.ext, where: 

n indicates the scale Q 2 to which data refers: n — i refers to the initial scale, the higher 
value of n = f refers to the final scale; if n is a number, then it refers to the n-th 
quark (ordered by increasing mass) production thresholds; 

i is the identificative of the distribution, reported in the third column of the table in 
subsection 2.5.4; 



ext is an extension chosen by the user. 



2.6 Running the code 

In this section we show some pdf plots obtained by running the code in many different 
situations. 

In Figg. (2.1 - 2.3), the evolution of up, down and gluon distributions at Q — 10 GeV 
and Q = 100 GeV in the unpolarized case is plotted at LO, NLO and NNLO. As initial 
distributions we have used both MRST and Alekhin distributions at the scale Qq = 
1.25 GeV 2 . Some details about the inputs are shown in the Table (2.1). 
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(f) Alekhin NNLO 



Figure 2.1: Evolution of the quark up distribution xu versus x at different Q values. All 
perturbative orders and both input models are shown. 
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Figure 2.2: Evolution of the quark down distribution xd versus x at different Q values. 
All perturbative orders and both input models are shown. 
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Figure 2.3: Evolution of the gluon distribution xg versus x at different Q values. All 
perturbative orders and both input models are shown. 
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Chapter 3 



NNLO precision studies of the Higgs 
total cross section 

3.1 Introduction 

One of the main objective at the LHC will be the search of the Higgs particle, which is 
responsible for the mechanism of mass generations for the fermions and the gauge bosons 
within the Standard Model. At the moment one of the most widely accepted ways to 
generate a mass for these particles relies on a scalar sector with a single Higgs field and a 
symmetry breaking potential. While the Lagrangian is invariant under a local symmetry 
of SU(2) x U(l) and this, of course, requires the same symmetry for the scalar sector 
as well, the vacuum of the theory, which is a stable state of minimum energy, doesn't 
share the same symmetry. The presence of a minimal Higgs sector is also specific of the 
Standard Model, though many extensions of the Standard Model, supersymmetric and 
not, require a far more complex Higgs structure to give mass to both up-type and down- 
type quarks and to the leptons by some trilinear interactions termed Yukawa couplings. 
Our objective, in this chapter, is to present very accurate predictions for the total Higgs 
cross section by combining results available in the literature, specifically the computation 
of NNLO hard scatterings for the production of a scalar and a pseudoscalar Higgs in pp 
collisions. The analytical expressions are available through the work of Ravindran, Smith 
and van Neerven [107]. These have been combined with the NNLO evolution of the parton 
distributions in order to generate very accurate cross sections for the Higgs search at the 
LHC. 

After the submission of this thesis, the preliminar work presented in this chapter has 
been completed and published in [31]. 
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3.2 Higgs detection at hadron colliders 

We summarize some generalities on the Standard Model and the mechanism of mass 
generation, with an emphasis on their phenomenological implications to hadron colliders. 

The electroweak sector of the standard model is a chiral theory SU{2)l x U(1)y gauge 
theory containing 3 SU(2) gauge bosons, W^, and one U(l) gauge boson, B^. Left-handed 
and right-handed quarks couple differently to the 3 gauge bosons of the weak isospin and 
to B^. The Lagrangian is given by 

C KE = -- A W; V W^ - -B^ (3.1) 

where 



B^v = dyB^ — d^By . (3.2) 
A complex scalar field, doublet of SU (2) is introduced and defined as 

• = -L(* + i M , (3.3) 

and with a scalar potential given by 

V{$) = fi 2 | | +A^| |^) 2 , (3.4) 

(A > 0). This turns out to be the most general renormalizable and 577(2) invariant 
potential allowed. The minimum of the potential is now a valley of equivalent vacua, 
among which we choose a vacuum by a suitable parameterization. One possibility is 
given by 

<*> = 4 f ° ! (3-) 



V2 



v 



with a U(1)y hypercharge Y$ = 1. The electromagnetic charge is Q = T 3 + ^ and once 
can easily check that this state has not electric charge since 

Q{&) = (3.6) 



The vacuum expectation value of Eq. 3.5 allows to leave the U(l) em intact by making one 
linear combination of the third component of 577(2) and U(l)y gauge symmetric, with Q 
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identified with the generator of U(l) em , according to the braking pattern 

SU(2) L x U(l)y — > U(1)em- The masses of the gauge bosons are generated by the 
kinetic part of the Higgs Lagrangian 

C s = (D^(D^) - (3.7) 
where r, are the Pauli matrices. The explicit form of the covariant derivatives is 

D (i = d ll + i^T-W lt + i^B li . (3.8) 
In unitary gauge the scalar field can be written as 

*=^=( ° ) (3-9) 
y/2 \ v + H J V ' 

and the contribution to the gauge boson masses from the scalar kinetic energy term are 

2 



5(0 v^gr-W. + fa) ° . (3.10) 



To summarize, the masses of the massive gauge bosons 

1 



3 



-g'B, + gWl 



Vg 2 + g' 2 
v g + g 



are given by 



H4-2 ^22 

M w = 13 v 



Ml = \(g 2 + g' 2 y 

M 1 = 0. (3.12) 
and the coupling constants satisfy relations 

e = g sin 6 W 

e = g'cos9 w (3.13) 



with thetaw being the Weinberg angle. Notice that in the unitary gauge a counting of 
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the degrees of freedom for the Higgs scalar and the gauge bosons remains the same, as 
expected, without the appearance of unphysical Goldstone bosons. Alternative parame- 
terizations, a la Kibble 




$ = -/H , „ ■ ( 3 - 14 ) 



let the Goldstone modes appear from the scalar sector. The appearance of this modes 
is a gauge variant effect. Non unitary gauges are important especially in the context of 
testing the unitarity of the model. In the Standard Model, there are three Goldstone 
bosons, Q = (w ± ,z), with masses Mw and Mz in the Feynman gauge. 

In addition to giving the W and Z bosons their masses, the Higgs boson is also 
responsible for giving the mass to the fermions. Thegauge invariant Yukawa coupling of 
the Higgs boson to fermions is given by trilinear terms of the form 

C f = -\ d Q L <S>d R + h.c. , (3.15) 

where the left handed SU(2) fermion doublet is constructed by acting with the chiral 
Lorenz projectors on both components 




Ql = , • (3.16) 



This gives the effective coupling 

^d-^=(u L , d L ) f ° \d R + h.c. (3.17) 
V 2 \ v + H 

which can be seen to yield a mass term for the down quark if we make the identification 

A, = ^. (3.18) 

V 

The mass term for the up quark is obtained by giving a value expectation value to the field 
$ c = — ?r 2 $* which is an 577(2) doublet and allows to write down the 577(2) invariant 
coupling 

\ u Q L <5> c u R + h.c. (3.19) 

which is the mass term for the up quark. For the upper components of the lepton SU (2) 
doublets the procedure is the same, though the generation of a mass term for the neutrino 
is still a problem within the Standard Model. The neutrino appears only in its left- 
component and therefore has no Yukawa coupling. 
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t" > h 



Figure 3.1: The leading order diagram for Higgs production by gluon fusion 

For the multi- family case, the Yukawa couplings, and A u , become matrices which 
are valued in family space (Np x Np) with Np being the number of families. If we 
introduce a family mixing matrix U and define g<i and gy to be the SU(2) and U(l)y 
coupling constant the final Lagrangian describing the interactions of the leptons with the 
massive gauge bosons is 



a *- ' :> i ; < ,-■ — IflV 7^ fi'A 7^7 ) e 

2 costly v 7 zcostV v ' 

(3.20) 

where 

= TiS - 2 sin 2 9 f A - Z ° = (3.21) 



while for the quarks one obtains 
C q = -e (^ua'Ui - ^di7~di) A M 

2 cos 6V v 7 2 cos 6V v 7 

(3.22) 

A complete set of Feynman rules describing the interaction of the Higgs with the 
massive states of the Standard Model in both unitary and non unitary gauges, which 
involve the goldstones, can be found in [20]. The Higgs field, being responsible for the 
mechanism of mass generation, can be radiated off by any massive state and its coupling 
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Figure 3.2: A typical NLO diagram for Higgs production by gluon fusion 



is proportional to the mass of the same state. At the LHC one of the golden plated modes 
to search for the Higgs is its production via the mechanism of gluon fusion. The leading 
order contribution is shown in fig. (3.1) which shows the dependence of the amplitude 
through the quark loop. Most of the contribution comes from the top quark, since this 
is the heaviest and has a larger coupling to the Higgs field. NLO and NNLO corrections 
have been computed in the last few years by various groups. A typical NLO correction is 
shown in fig. (3.2). From the computation of the LO diagram one gets 

Ass - H) = -^(^M - pY .) ^M JImU ^^ 1 (3 ' 23) 

where a and b are color indices and p and q are the momenta of the two gluons. When 
the fermion in the loop is much heavier than the Higgs boson, m» M H the expression 
above becomes 

A(gg -> H) ^ m>>MlI -^(^y-^^M?)- ( 3 - 24 ) 

which apparently does not depend on the quark mass. This implies that if there are 
additional generators, then the contribution of this process is proportional to the number 
of additional generations, providing an important window into new physics beyond the 
Standard Model. The resonant production of a heavy Higgs is instead given by the formula 

7T 2 M 2 

H " ^ H) = Mdf [H ~* 99)6(1 " ~f ] - (3 - 25) 
relating the cross section for Higgs production via gluon fusion to the corresponding decay 
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with 



rate F(H — > gg). At parton level we obtain 

a(gg - H) = | l(^f) | 2 5{s - M%) (3.26) 

where v^§ is the energy in the gluon -gluon center of mass and the integral I is defined by 

3.2.1 The effective Lagrangian 

An important feature of the result for Higgs boson production from gluon fusion is that it 
is independent of the heavy quark mass for a light Higgs boson. Eq. 3.24 can be derived 
from an effective vertex 

PF = id (3 - 28) 

being the contribution of heavy fermion loops to the QCD beta function. The effective 
Lagrangian can be used to compute the radiative corrections in the gluon sector. The 
correction in principle involve 2-loop diagrams. However, using the effective vertices from 
Eq. 3.28, the O(a^) corrections can be found from a 1-loop calculation. In this way the 
contribution from the top quark loop is shrunk to a point. A discussion of the NNLO 
approach to the computation of the gluon fusion contributions to Higgs production has 
been presented in [107], work to which we refer for more details. We recall that in this 
work the authors present a study for both scalar and pseudoscalar Higgs production, the 
pseudoscalar appearing in 2-Higgs doublets models. In these extended models there is 
a mixing between the two scalars which is usually parameterized in terms of the ratio 
between the two vacuum expectation values at the minimum of the potential i>i/i>2 = 
tan/3. Diagonalization of the mass matrix for the Higgs at the minimum introduces scalars 
and pseudoscalars interactions between the various Higgs and the quarks, as shown from 
the structure of the operator Oi- Here we will briefly summarize their results. 

In the large top-quark mass limit the Feynman rules for scalar Higgs production (H) 
can be derived from the effective Lagrangian density 

Cf ff = G u ^(x)0(x) with 0(x) = -^G%(x)G a ^(x), (3.29) 
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whereas the production of a pseudo-scalar Higgs (A) is obtained from 



-eff 



G A 1 (x)+G A 2 (x) 



with 



1 



1 (x) = ~e^G^G^(x), 



2 (x) = -- Qi( x ) 7m 75 Qi(x) , 

A i=i 



(3.30) 



where $ H (x) and $ A (a;) represent the scalar and pseudo-scalar fields respectively and nf 
denotes the number of light flavours. G%" is the field strength of QCD and the quark field 
are denoted by With an appropriate normalization the constants turn out to be given 

by 

G B = -2^a s {^)Gfr B F B {r B )C B ^a s {^),^ , 



G A 



and a s (f/%) is defined by 



a s (fi 2 r )C F Q -3 ln£| 



Ga , 



(3.31) 



2, '>,(//;) 



ill 



(3.32) 



where (^(/i 2 .) is the running coupling constant and \i r denotes the renormalization scale. 
Gf is the Fermi constant and the functions F B are given by 

F H (r) = 1 + (1 - r) /(r) , F A (r) = /(r) cot /? , 
4 m? 

/(r) = arcsin 2 — p , for r > 1 , 



/00 = -7 h 



1 , 1- 



4 I l + v^ 



+ 7r i for r < 1 , 



(3.33) 



where cot (5 denotes the mixing angle in the Two-Higgs-Doublet Model. Further m and 
m t denote the masses of the (pseudo-) scalar Higgs boson and the top quark respectively. 
The coefficient functions C B originate from the corrections to the top-quark triangular 
graph provided one takes the limit m t — > oo. The coefficient functions are known up to 
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order a 2 s and are given by 



M^)A)=l + ai 5) (/4) 



mi 



5 Ca — 3Cp 



2 F 



m CACF+ ]^ cl .i CFT! Ji CATf+ i 7C 



-IICaCf) \n^ + n f T f ( -5C F - ^- C A + SC F In J 



C A (W), A<; 



1, 



(3.34) 
(3.35) 



where is given in the five-flavour number scheme. 

Using the effective Lagrangian approach one can calculate the total cross section of 
the reaction 



H 1 (P 1 )+H 2 (P 2 )^B+' X', 



(3.36) 



where Hi and H 2 denote the incoming hadrons and X represents an inclusive hadronic 
state. The total cross section is given by 



Ctot 



( N2 ~ 1) a,b= q , q , 9 J * J ^ 



I X 77T \ 

xA a6iB , B = H, A 

\x l x 2 fx 2 I 



m 
~S~ 



with x = , S = (Pi + P 2 f , p 2 5 = m 2 



(3.37) 



where the factor l/(iV 2 — 1) is due to the average over colour and the parton distributions 
f a (y,fi 2 ) (a,b = q,q,g) depend on the mass factorization/renormalization scale \i. A^b 
denotes the partonic hard scattering coefficient computed with NNLO accuracy. 



3.3 Results 

The use of our results on the NNLO evolution of the parton distributions together with 
the results of [107] allows us to provide accurate predictions for the total cross section 
for Higgs production. Here we summarize our numerical results which are illustrated in 
Figs. (3.3 - 3.6) and in Tables (3.1 - 3.4). The sets of distributions compared are those of 
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MRST and Alekhin, currently the only ones available at NNLO. The plots refer to center 
of mass energies which are typical at the LHC, 14 TeV being the largest achievable in 
a not so distant future. The pseudoscalar cross section is systematically larger than the 
scalar one and both drop as the mass of the Higgs particles increases. The interval of 
variability which has been considered for this parameter is light-to-heavy (100 GeV to 300 
GeV). The factorization scale and the renormalization scales have been chosen to coincide. 
For a given value of the Higgs mass, the cross section raises considerably with energy by 
a factor approximately 50-70 as we move from the lower toward the higher value of S, 
the center of mass energy. Other distinct features of this study is the dependence of the 
result on the initial model chosen for the evolution. Finally the inclusion of the NNLO 
corrections, in all models, causes an increase in the actual values of the cross sections 
compared to the NLO results. A possible extension of this study is the inclusion of 
renormalization effects in the evolution and in the hard scatterings, which can be studied 
along the lines discussed in chapter 2. 




100 150 200 250 300 100 150 200 250 300 

m H [GeV] m H [GeV] 

(a) MRST (b) Alekhin 



Figure 3.3: Total cross section as a function of Higgs mass at \/S = 14 TeV in pp collision, 
for a scalar Higgs boson. 



Chapter 3. NNLO precision studies of the Higgs total cross section 



83 





LO 

NLO — - 
MNLO 











200 
m H [GeV] 







( ) 

NLO — 
MNLO 















200 
m H [GeV] 



(a) MRST 



(b) Alekhin 



Figure 3.4: Total cross section as a function of Higgs mass at = 2 TeV in pp collision, 
for a scalar Higgs boson. 




LO 
NLO 
NNLO 



150 200 

m H [GeV] 




(a) MRST 



(b) Alekhin 



Figure 3.5: Total cross section as a function of Higgs mass at = 14 TeV in pp collision, 
for a pseudoscalar Higgs boson. 
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Figure 3.6: Total cross section as a function of Higgs mass at \fS = 2 TeV in pp collision, 
for a pseudoscalar Higgs boson. 





LO 


NLO 


NNLO 


m H 


MRST 


Alekhin 


MRST 


Alekhin 


MRST 


Alekhin 


100 


33.12 


32.02 


53.87 


54.89 


57.31 


70.46 


110 


27.86 


26.97 


45.99 


46.69 


49.24 


59.81 


120 


23.75 


23.01 


39.74 


40.21 


42.81 


51.42 


130 


20.48 


19.86 


34.69 


35 


37.58 


44.69 


140 


17.83 


17.31 


30.57 


30.77 


33.29 


39.22 


150 


15.67 


15.23 


27.17 


27.28 


29.74 


34.73 


160 


13.89 


13.5 


24.33 


24.37 


26.75 


30.99 


170 


12.39 


12.06 


21.92 


21.92 


24.22 


27.84 


180 


11.13 


10.84 


19.87 


19.84 


22.05 


25.16 


190 


10.06 


9.8 


18.12 


18.05 


20.19 


22.88 


200 


9.138 


8.912 


16.61 


16.52 


18.58 


20.91 


210 


8.35 


8.149 


15.3 


15.19 


17.18 


19.22 


220 


7.673 


7.492 


14.17 


14.05 


15.97 


17.75 


230 


7.089 


6.925 


13.18 


13.06 


14.92 


16.49 


240 


6.584 


6.435 


12.33 


12.2 


14.01 


15.39 


250 


6.148 


6.012 


11.6 


11.46 


13.22 


14.44 


260 


5.771 


5.646 


10.96 


10.81 


12.53 


13.62 


270 


5.446 


5.33 


10.41 


10.26 


11.94 


12.91 


280 


5.169 


5.061 


9.94 


9.786 


11.44 


12.31 


290 


4.935 


4.834 


9.548 


9.39 


11.03 


11.81 


300 


4.743 


4.648 


9.23 


9.069 


10.7 


11.4 



Table 3.1: Total cross section as a function of Higgs mass at yfS = 14 TeV in pp collision, 
for a scalar Higgs boson. The mass of the Higgs boson is in GeV and the cross sections 
are in pb. 
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LO 


NLO 


NNLO 


m H 


MRST 


Alekhin 


MRST 


Alekhin 


MRST 


Alekhin 


100 


0.6434 


0.6285 


1.528 


1.47 


1.974 


1.93 


110 


0.4838 


0.4737 


1.172 


1.126 


1.534 


1.474 


120 


0.37 


0.3632 


0.9136 


0.8773 


1.211 


1.145 


130 


0.2871 


0.2826 


0.7219 


0.6931 


0.9682 


0.9021 


140 


0.2256 


0.2228 


0.5775 


0.5545 


0.783 


0.7196 


150 


0.1793 


0.1777 


0.467 


0.4486 


0.6399 


0.5806 


160 


0.1439 


0.1432 


0.3812 


0.3665 


0.5277 


0.4731 


170 


0.1165 


0.1164 


0.3139 


0.302 


0.4387 


0.3888 


180 


0.09506 


0.09537 


0.2603 


0.2507 


0.3673 


0.3219 


190 


0.07814 


0.07874 


0.2174 


0.2095 


0.3095 


0.2685 


200 


0.06466 


0.06546 


0.1828 


0.1763 


0.2625 


0.2254 


210 


0.05385 


0.05478 


0.1546 


0.1493 


0.2239 


0.1904 


220 


0.04514 


0.04615 


0.1316 


0.1272 


0.1922 


0.1618 


230 


0.03806 


0.03912 


0.1127 


0.109 


0.1659 


0.1384 


240 


0.03228 


0.03336 


0.09704 


0.09395 


0.144 


0.119 


250 


0.02754 


0.02862 


0.08408 


0.08146 


0.1256 


0.1029 


260 


0.02363 


0.0247 


0.07326 


0.07102 


0.1102 


0.08954 


270 


0.0204 


0.02144 


0.0642 


0.06227 


0.09728 


0.07833 


280 


0.01771 


0.01872 


0.0566 


0.05492 


0.08633 


0.06893 


290 


0.01547 


0.01646 


0.05022 


0.04874 


0.07709 


0.06103 


300 


0.01361 


0.01457 


0.04486 


0.04355 


0.0693 


0.05441 



Table 3.2: Total cross section as a function of Higgs mass at V~S = 2TeV in pp collision, 
for a scalar Higgs boson. The mass of the Higgs boson is in GeV and the cross sections 
are in pb. 
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LO 


NLO 


NNLO 


m H 


MRST 


Alekhin 


MRST 


Alekhin 


MRST 


Alekhin 


100 


75.83 


73.31 


124.3 


126.7 


132.6 


163.1 


110 


64.04 


61.98 


106.5 


108.1 


114.3 


139 


120 


54.82 


53.12 


92.44 


93.54 


99.84 


120 


130 


47.5 


46.06 


81.09 


81.81 


88.08 


104.8 


140 


41.59 


40.37 


71.85 


72.3 


78.46 


92.47 


150 


36.77 


35.72 


64.23 


64.49 


70.5 


82.37 


160 


32.79 


31.89 


57.88 


58 


63.84 


73.98 


170 


29.47 


28.68 


52.54 


52.54 


58.22 


66.95 


180 


26.68 


25.99 


48.01 


47.93 


53.43 


60.99 


190 


24.32 


23.7 


44.15 


44 


49.35 


55.94 


200 


22.32 


21.77 


40.87 


40.65 


45.86 


51.64 


210 


20.61 


20.12 


38.05 


37.79 


42.87 


47.96 


220 


19.17 


18.72 


35.66 


35.36 


40.33 


44.84 


230 


17.95 


17.53 


33.63 


33.31 


38.18 


42.2 


240 


16.92 


16.54 


31.94 


31.59 


36.39 


39.99 


250 


16.07 


15.71 


30.54 


30.16 


34.91 


38.15 


260 


15.37 


15.04 


29.41 


29.02 


33.74 


36.68 


270 


14.83 


14.51 


28.55 


28.14 


32.87 


35.55 


280 


14.43 


14.13 


27.97 


27.53 


32.3 


34.76 


290 


14.2 


13.91 


27.68 


27.22 


32.07 


34.35 


300 


14.15 


13.87 


27.74 


27.26 


32.25 


34.38 



Table 3.3: Total cross section as a function of Higgs mass at v^S* = 14TeV in pp collision, 
for a pseudoscalar Higgs boson. The mass of the Higgs boson is in GeV and the cross 
sections are in pb. 
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LO 


NLO 


NNLO 


m H 


MRST 


Alekhin 


MRST 


Alekhin 


MRST 


Alekhin 


100 


1.473 


1.439 


3.526 


3.393 


4.579 


4.477 


110 


1.112 


1.089 


2.715 


2.609 


3.572 


3.432 


120 


0.8541 


0.8385 


2.125 


2.041 


2.831 


2.677 


130 


0.6659 


0.6556 


1.688 


1.62 


2.274 


2.119 


140 


0.5261 


0.5196 


1.357 


1.303 


1.849 


1.699 


150 


0.4205 


0.4168 


1.104 


1.06 


1.52 


1.379 


160 


0.3397 


0.338 


0.9071 


0.8719 


1.262 


1.131 


170 


0.2771 


0.2768 


0.7522 


0.7236 


1.056 


0.9361 


180 


0.2279 


0.2287 


0.6288 


0.6055 


0.8914 


0.7813 


190 


0.189 


0.1904 


0.5298 


0.5106 


0.7578 


0.6572 


200 


0.1579 


0.1599 


0.4497 


0.4338 


0.6489 


0.557 


210 


0.1329 


0.1352 


0.3844 


0.3712 


0.5595 


0.4756 


220 


0.1128 


0.1153 


0.3311 


0.32 


0.4859 


0.4091 


230 


0.09636 


0.09904 


0.2873 


0.2779 


0.4249 


0.3544 


240 


0.08297 


0.08573 


0.2512 


0.2432 


0.3744 


0.3094 


250 


0.07198 


0.07479 


0.2213 


0.2144 


0.3322 


0.2722 


260 


0.06295 


0.06578 


0.1965 


0.1905 


0.2971 


0.2413 


270 


0.05553 


0.05837 


0.176 


0.1708 


0.268 


0.2157 


280 


0.04945 


0.05228 


0.1592 


0.1545 


0.2439 


0.1947 


290 


0.04452 


0.04735 


0.1455 


0.1412 


0.2244 


0.1776 


300 


0.0406 


0.04346 


0.1348 


0.1308 


0.2091 


0.1642 



Table 3.4: Total cross section as a function of Higgs mass at V~S = 2 TeV in pp collision, 
for a pseudoscalar Higgs boson. The mass of the Higgs boson is in GeV and the cross 
sections are in pb. 
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Chapter 4 



The kinetic interpretation of the 
DGLAP Equation, its Kramers- Moyal 
expansion and positivity of helicity 
distributions 



4.1 Background on the topic 



In this chapter and in the chapter that follows we elaborate on what has come to be 
known as the "kinetic description" of the dynamics of partons under the action of the 
renormalization group. In the DGLAP evolution the parton densities tend to be supported 
by smaller and smaller values of the Bjorken variable x as the energy variable log(Q) 
increases. Building on previous observations by Teryaev and by Collins and Qiu we 
clarify that the DGLAP evolution can be interpreted as a kinetic master equation which 
can be expanded formally to generate a hierarchy of new equations using the Kramers- 
Moyal approach. Arrested to the second order, this expansion generates a Fokker-Planck 
approximation to the DGLAP evolution. The issue of positivity of the kernels is important 
in order for this picture to hold. The LO kernels have this property, while at NLO a formal 
proof of positivity is more intricate and requires a numerical analysis to be supported. 
The hierarchy of equations generated by the Kramers-Moyal expansion could be studied 
independently, possibly numerically in order to gain more insight into the property of this 
expansion. One of the consequence of the kinetic interpretation, as shown in this chapter 
is the proof is that parton dynamics is reminscent of a "downward biased random walk". 
These issues are analyzed in detail. 
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4.2 Introduction 

According to a rederivation - due to Collins and Qiu - the DGLAP equation can be rein- 
terpreted (in leading order) in a probabilistic way. This form of the equation has been 
used indirectly to prove the bound \ Af(x,Q)\ < f(x,Q) between polarized and unpolar- 
ized distributions, or positivity of the helicity distributions, for any Q. We reanalyze this 
issue by performing a detailed numerical study of the positivity bounds of the helicity 
distributions. We also elaborate on some of the formal properties of the Collins-Qiu form 
and comment on the underlying regularization, introduce a Kramers-Moyal expansion of 
the equation and briefly analyze its Fokker-Planck approximation. These follow quite nat- 
urally once the master version is given. We illustrate this expansion both for the valence 
quark distribution qy and for the transverse spin distribution h\ and analyze the role of 
one well known inequality from this perspective. In fact, an interesting constraint relating 
longitudinally polarized, unpolarized and transversely polarized distributions is Soffer's 
inequality, which deserves a special attention, since has to be respected by the evolution 
to any order in a s . Some tests of the inequality have been performed in the near past, 
bringing support to it. However, other inequalities are supposed to hold as well. 

In this work we perform a NLO analysis of an inequality which relates longitudinally 
polarized distributions and unpolarized ones. The inequality can be summarized in the 
statement that helicity distributions (positive and negative) for quarks and gluons have 
to be positive. The inequality states that 



where the ± refers to the the possible values of the helicities of quarks and gluons. The 
statement is supposed to hold, at least in leading order, for any Q. To analyze the renor- 
malization group evolution of this relation, especially to next-to-leading order, requires 
some effort since this study involves a combined study of the (longitudinally) polarized 
and unpolarized evolutions. In this work we present a complete NLO study of the evo- 
lution equations starting directly from the helicity basis. Helicities are in fact the basic 
parton distributions from which other distributions can be built. 

Compared to other implementations, in our work we perform a NLO test of the pos- 
itivity of the helicity distributions using an ansatz due to Rossi [109] which reduces the 
evolution equations to an infinite set of recursion relations for some scale invariant coef- 
ficients. 

Various arguments to validate eq. (4.2) have been presented in the literature. From 




(4.1) 



or 



f ± (x,Q 2 )>0 



(4.2) 
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our perspective, an interesting one has been formulated by Teryaev and Collaborators 
who have tried to establish a link, to leading order, between evolution equations and their 
probabilistic interpretation in order to prove Soffer's inequality. Similar arguments hold 
also in the analysis of eq. (4.2). 

We should remark that a complete probabilistic picture exists only for the leading 
order unpolarized evolution [38] and the arguments of [21] are inspired by the fact that 
the subtraction terms (the x — 1 contributions in the expressions of the kernels, where 
x is Bjorken's variable), being positive, once they are combined with the bulk (x < 1) 
contributions give a form of the evolution equations which are diagonal in parton type 
and resemble "kinetic" equations. Our arguments, on this issue, are just a refinement of 
this previous and influential analysis. 

In the recent literature there has been some attention to this feature of the DGLAP 
evolution, limited to the nonsinglet sector, in connection with kinetic theory and the 
"dynamical renormalization group", in the words of ref. [23]. 

All the arguments, so far, go back to some important older work of Collins and Qiu who 
provided an interesting derivation of the (unpolarized) DGLAP equation using Mueller's 
formalism of cut diagrams. In their paper [38] the authors reinterpreted the DGLAP 
equation as a kinetic probabilistic equation of Boltzmann type. The authors gave no detail 
on some of the issues concerning the regularization of their diagrammatic expansion, on 
which we will elaborate since we need it for our accurate numerical analysis. In our work 
the Collins-Qiu form of the DGLAP equation is interpreted simply as a master equation 
rather than a Boltzmann equation, given the absence of a 2-to-2 scattering cross section 
in the probabilistic partonic interpretation. A master equation is governed by transition 
probabilities and various formal approximations find their way once this conceptual step is 
made. We illustrate, in the spirit of a stochastic approach to the DGLAP dynamics, how to 
extract standard differential equations of Kramers- Moyal type for the simplest nonsinglet 
evolutions, those involving valence distributions and transverse spin distributions. 

We show that the DGLAP dynamics can be described, at least in a formal way, by a 
differential equation of arbitrarily high order. Truncations of this expansion to the first 
few orders provide the usual link with the Fokker-Planck approximation, the Langevin 
equation and its path integral version 1 . The picture one should have in mind, at least in 
this approximation, is that of a stochastic (Brownian) dynamics of Bjorken's variable x in 
a fictitious time \og(Q), describing the evolution under the renormalization group (RG). 
In this interpretation the probability function is the parton distribution itself. 

This chapter is based on the paper [26]. 



1 For an example of this interplay between differential and stochastic descriptions we refer the reader 
to [12] 
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4.3 Master Equations and Positivity 

Let's start considering a generic 1-D master equation for transition probabilities w(x\x') 
which we interpret as the probability of making a transition to a point x given a starting 
point x' for a given physical system. The picture we have in mind is that of a gas of 
particles making collisions in 1-D and entering the interval (x, x + dx) with a probability 
w(x\x') per single transition, or leaving it with a transition probability w(x'\x). In general 
one writes down a master equation 

r ) = J dx' (w(x\x')f(x', t) - w(x'\x)f(x, r)) dx'. (4.3) 

describing the time r evolution of the density of the gas undergoing collisions or the motion 
of a many replicas of walkers of density f(x,r) jumping with a pre-assigned probability, 
according to taste. 

The result of Collins a Qiu, who were after a derivation of the DGLAP equation 
that could include automatically also the "edge point" contributions (or x=l terms of the 
DGLAP kernels) is in pointing out the existence of a probabilistic picture of the DGLAP 
dynamics. These edge point terms had been always introduced in the past only by hand 
and serve to enforce the baryon number sum rule and the momentum sum rule as Q, the 
momentum scale, varies. 

The kinetic interpretation was used in [21] to provide an alternative proof of Soffer's 
inequality. We recall that this inequality 

\hi(x)\ < q + (x) (4.4) 

famous by now, sets a bound on the transverse spin distribution hi(x) in terms of the 
components of the positive helicity component of the quarks, for a given flavour. The 
inequality has to be respected by the evolution. We recall that hi, also denoted by the 
symbol 

A T q(x, Q 2 ) = q1(x, Q 2 ) - q^x, Q 2 ), (4.5) 

has the property of being purely nonsinglet and of appearing at leading twist. It is 
identifiable in transversely polarized hadron-hadron collisions and not in Deep Inelastic 
Scattering (DIS), where can appear only through an insertion of the electron mass in the 
unitarity graph of DIS. 

The connection between the Collins-Qiu form of the DGLAP equation and the master 
equation is established as follows. The DGLAP equation, in its original formulation is 
generically written as 
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x y 1 



Figure 4.1: The constrained random walk of the parton densities 



dq(x,Q 2 ) f 1 dy 2 

where we are assuming a scalar form of the equation, such as in the nonsinglet sector. 
The generalization to the singlet sector of the arguments given below is, of course, quite 
straightforward. To arrive at a probabilistic picture of the equation we start reinterpreting 
r = log(Q 2 ) as a time variable, while the parton density q(x,r) lives in a one dimensional 
(Bjorken) x space. 

We recall that the kernels are defined as "plus" distributions. Conservation of baryon 
number, for instance, is enforced by the addition of edge-point contributions proportional 
to 5(1 — x). 

We start with the following form of the kernel 

P(z) = P(z) - 5(1 - z) f 1 P(z) dz, (4.7) 

Jo 

where we have separated the edge point contributions from the rest of the kernel, here 
called P(z). This manipulation is understood in all the equations that follow. The 
equation is rewritten in the following form 



d , v f 1 , a fx\ q(y,r) y* dy 6 ( y\ g(x,T ; 
-q(x, r) = I dyP (- j — J q -P —— (4.8) 

Now, if we define 



w(x\y) = ^-P(x/y) 6 -^^ (4.9) 
7T y 

(4.8) becomes a master equation for the probability function q(x,r) 

There are some interesting features of this special master equation. Differently from other 
master equations, where transitions are allowed from a given x both toward y > x and 
y < x, in this case, transitions toward x take place only from values y > x and leave the 
momentum cell (x,x + dx) only toward smaller y values (see Fig. (4.1). 

Clearly, this sets a direction of the kinetic evolution of the densities from large x values 
toward smaller-x values as r, the fictitious "time" variable, increases. 
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4.4. Probabilistic Kernels 



Probably this is the simplest illustration of the fact that parton densities, at large 
final evolution scales, are dominated by their small-x behaviour. As the "randomly moving 
partons" reach the x ~ region of momentum space, they find no space where to go, while 
other partons tend to pile up toward the same region from above. This is the picture of a 
random walk biased to move downward (toward the small-x region) and is illustrated in 
Fig. (4.1). 



4.4 Probabilistic Kernels 

We briefly discuss some salient features of the structure of the kernels in this approach 
and comment on the type of regularization involved in order to define them appropriately. 

We recall that unpolarized and polarized kernels, in leading order, are given by 

PPs = ^ = C,(^-l-* + ?,<l-*)) 
P(J) = 27) (* 2 + (1 - x) 2 )) 

p( o) = g i+a-*) 2 



where 



and 



99 



X 



P<® = 2N c ( Tr ^-^ + l-2 + x(l-x))+^-5(l-x) (4.11) 



; 1 — x)+ x 



Cf = ^4, T f = T R n f = \n h ft = ^-N c - |z> (4.12) 



APW = 2T f {2x-l) 
APW = C F (2-s) 

APW = 2iV c f— -2, + l)+5(l-4, (4-13) 



1 — x 



+ 



while the LO transverse kernels are given by 



A ^«= c 4a^)- + - 2+ l S{1 - x) )- <4 ' 14) 
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The unpolarized kernels should be compared with the Collins-Qiu form 

Pqq = Iqq ~ S(l ~ x) ( dz^ qq 

JO 

P gg = 7<?<? ~ ( n / J dz lig + 2 J dzi ") 5 ^ ~~ ^ 



Pqg 7gg 
Pgq = ^gi 



(4.15) 



where 



Iqq = C F {—-l-X 

Iqg = ( 2x - 1) 

Igq = C F (2-X) 

laa = 2N c (-^— + --2 + x{l-x) 
V 1 - x x 



(4.16) 



These kernels need a suitable regularization to be well defined. Below we will analyze the 
implicit regularization underlying eq. (4.15). One observation is however almost immedi- 
ate: the component P gg is not of the form given by eq. (4.7). In general, therefore, in the 
singlet case, the generalization of eq. (4.7) is given by 

P(x) = PAx) - 5(1 - x) f 1 P 2 {z)dz (4.17) 

JO 

and a probabilistic interpretation is more complex compared to the nonsinglet case and 
has been discussed in the original literature [38]. 



4.5 Convolutions and Master Form of the Singlet 

Distributions are folded with the kernels and the result rearranged in order to simplify the 
structure of the equations. Since in the previous literature this is done in a rather involute 
way [53] we provide here a simplification, from which the equivalence of the various forms 
of the kernel, in the various regularizations adopted, will be apparent. All we need is the 
simple relation 



I 



dy f, 1 \ f 1 dy yf{y)-xf{x) 

f(x y)= + log(l - x)f(x) (4.18) 

Jx v v — x 



* y(i -y)+ J* y y-x 
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4.5. Convolutions and Master Form of the Singlet 



in which, on the right hand side, regularity of both the first and the second term is explicit. 
For instance, the evolution equations become 

dq nn f dy yq(y) - xq(x) f 1 dy 3 

di^m = 2Cf Jj y -x +2Cf log(1 - x) q{x) -Lj {1+z) q{y) + 2 CFq{x) 

+n f f ^ (z 2 + (1 - zf) g(y) 
dq „ f 1 dy 1 + (1 - z) 2 , . „ T f 1 dy y f (y) - xf (x) , s 

+2JV c log(l - x)g(x) + 2N C f ^ (- - 2 + *(1 - z)) j(t/) + ^g(x) 

Jx y \z / I 



(4.19) 



with z = x/y. The same simplified form is obtained from the probabilistic version, having 
defined a suitable regularization of the edge point singularities in the integrals over the 
components •jff in eq. (4.16). The canonical expressions of the kernels (4.13), expressed in 
terms of "+" distributions, can also be rearranged to look like their equivalent probabilistic 
form by isolating the edge-point contributions hidden in their "+" distributions. We get 
the expressions 



and 

AP™ = 2N C - 2,+ l) - (2^/; ^ - |) ,<1 - *), (4.2!) 

the other expressions remaining invariant. In appendix A we provide some technical 
details on the equivalence between the convolutions obtained using these kernels with the 
standard ones. 

A master form of the singlet (unpolarized) equation is obtained by a straightforward 
change of variable in the decreasing terms. We obtain 

^: = j x < ^-i qq (x/y)q(y) - J q < ^-i qq (y/x)q(x) 
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% = I ^-lgg^ly)-nf f l qg {y/x)g{x) 
dr Jx y Jo 

"9 i l 99 {y/x)g{x) + / -f>y gq (x/y)q(y) (4.22) 

A J A Jx y 

with a suitable (unique) cutoff A needed to cast the equation in the form (4.19). A 
discussion of this aspect is left in appendix B. The (regulated) transition probabilities are 
then given by 

w qq {x\y) = l qq {x/y) 9 y> x 9 y <l — ^_ 

w qq {y\x) = -f qq (y/x) — < X ^)^ >Q ) 

9(y > x)9(y < 1 - A) 



w gg (x\y) = i gg {x/y)- 



y 



W qq {y\x) = ("./-,/„!//.'•) -~gg(<J-<-^ 



2 ' yyw ' V x 

Wg q {y\x) = Iggix/y) 9 ^ > X ^ < 1 — 

Wgq{x\y) = 0, 

(4.23) 

as one can easily deduct from the form of eq. (4.10). 



4.6 A Kramers- Moyal Expansion for the DGLAP Equa- 
tion 

Kramers-Moyal (KM) expansions of the master equations (backward or forward) are some- 
times useful in order to gain insight into the master equation itself, since they may provide 
a complementary view of the underlying dynamics. 

The expansion allows to get rid of the integral which characterizes the master equation, 
at the cost of introducing a differential operator of arbitrary order. For the approximation 
to be useful, one has to stop the expansion after the first few orders. In many cases this 
turns out to be possible. Examples of processes of this type are special Langevin processes 
and processes described by a Fokker-Planck operator. In these cases the probabilistic in- 
terpretation allows us to write down a fictitious Lagrangian, a corresponding path integral 
and solve for the propagators using the Feynman-Kac formula. For definitess we take the 
integral to cover all the real axis in the variable x' 
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4.6. A Kramers-Moyal Expansion for the DGLAP Equation 



— — q(x, t) = / dx' (w(x\x')q(x' i r) — w(x'\x)q(x, r)) dx'. (4.24) 

OT J — oo 

As we will see below, in the DGLAP case some modifications to the usual form of the 
KM expansion will appear. At this point we perform a KM expansion of the equation in 
the usual way. We make the substitutions in the master equation y — > x — y in the first 
term and y — > x + y in the second term 

d f°° 

—q(x,T)— / dy (w(x\x — y)q(x — y, r) — w(x + y\x)q(x, r)) , (4-25) 

OT J-oo 

identically equal to 

_g(x,r)=/ dy(w(x + y-y'\x-y')q(x-y',T) - w(x + y'\x)q(x,T)) , (4.26) 
or J-oo 

with y = y' . First and second term in the equation above differ by a shift (in —y') and 
can be related using a Taylor (or KM) expansion of the first term 



—q(x,r) = J_^dyY^ ~ } ^ (w(x + y\x)q(x, r)) (4.27) 

where the n = term has canceled between the first and the second contribution coming 
from (4.26). The result can be written in the form 

f) 00 (— v ) n B n 

9;fa r) - E ~^Td^ M*M X > ')) ( 4 - 28 ) 

where 

/oo 
dy(y-x) n w(y\x). (4.29) 
-oo 

In the DGLAP case we need to amend the former derivation, due to the presence of 
boundaries (0 < x < 1) in the Bjorken variable x. For simplicity we will focus on the 
nonsinglet case. We rewrite the master equation using the same change of variables used 
above 



d f 1 f x 

— q ( X)T ) = J dyw(x\y)q(y,r) - dyw(y\x)q(x,r) 

ra(x) r-x 

— / dywix + y\x) * q(x, r) + / dywix + y \x)q(x, r), (4.30) 

JO JO 

where we have introduced the simplest form of the Moyal product 2 

2 A note for noncommutative geometers: this simplified form is obtained for a dissipative dynamics 
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w(x + y\x) * q(x) = w(x + y\x)e V ( x+ x ^q(x,r) (4-31) 
and a(x) — x — 1. The expansion is of the form 

—q(x,T)= dyw(x + y\x)q(x,r) -^2 d V ^—d x n (w(x + y\x)q(x,T))(4.32) 

OT Ja(x) n=l"' n - 

which can be reduced to a differential equation of arbitrary order using simple manipula- 
tions. We recall that the Fokker-Planck approximation is obtained stopping the expansion 
at the second order 

d 1 

— q(x, t) = a (x) - d x (ai(x)q(x)) + -d 2 x (a 2 (x)q(x, r)) (4.33) 

with 

a n (x) = J dyy n w(x + y,x) (4.34) 

being moments of the transition probability function w. Given the boundary conditions 
on the Bjorken variable x, even in the Fokker-Planck approximation, the Fokker-Planck 
version of the DGLAP equation is slightly more involved than Eq. (4.33) and the coeffi- 
cients a n (x) need to be redefined. 



4.7 The Fokker-Planck Approximation 

The probabilistic interpretation of the DGLAP equation motivates us to investigate the 
role of the Fokker-Planck (FP) approximation to the equation and its possible practical 
use. We should start by saying a word of caution regarding this expansion. 

In the context of a random walk, an all-order derivative expansion of the master 
equation can be arrested to the first few terms either if the conditions of Pawula's theorem 
are satisfied -in which case the FP approximation turns out to be exact- or if the transition 
probabilities show an exponential decay above a certain distance allowed to the random 
walk. Since the DGLAP kernels show only an algebraic decay in x, and there isn't any 
explicit scale in the kernel themselves, the expansion is questionable. However, from a 
formal viewpoint, it is still allowed. With these caveats in mind we proceed to investigate 
the features of this expansion. 

We redefine 



when the p's of phase space are replaced by constants. Here we have only one variable: x 
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4.7. The Fokker-Planck Approximation 



r—x 

do(x) = / dyw(x + y\x)q(x, r) 

Ja(x) 
ra(x) 

a n (x) = / dyy n w(x + y\x)q(x,T) 
Jo 

rce(x) 

a n (x) = / dyy n d x n (w(x + y\x)q(x,r)) n = l,2, ... (4.35) 
Jo 



For the first two terms (n = 1, 2) one can easily work out the relations 

di(x) = d x ai(x) — a(x)d x a(x)w(x + a(x)\x)q(x, r) 

d 2 (x) = d x \a2(x) — 2a(x)(d x a(x)) 2 w(x + a(x)\x)q(x, r) — a(x) 2 d x a(x)d x (w(x + a(x)\x)q(x, r)) 
-a 2 (x)d x a(x)d x (w(x + y\x)q(x, r)) \ y=a (x) (4.36) 

Let's see what happens when we arrest the expansion (4.32) to the first 3 terms. The 
Fokker-Planck version of the equation is obtained by including in the approximation only 
d n with n — 0, 1, 2. 

The Fokker-Planck limit of the (nonsinglet) equation is then given by 

d 1 

— q(x, r) = a (x) + a x (x) - -a 2 (x) (4.37) 

which we rewrite explicitly as 

d . , _ /85 3 13 10 12 ol , . 

A similar approach can be followed also for other cases, for which a probabilistic 
picture (a derivation of Collins-Qiu type) has not been established yet, such as for hi. We 
describe briefly how to proceed in this case. 

First of all, we rewrite the evolution equation for the transversity in a suitable master 
form. This is possible since the subtraction terms can be written as integrals of a positive 
function. A possibility is to choose the transition probabilities 

wi[x\y] = ^ ( - 2 )9(y > x)9(y < 1) 

w 2 [y\x\ = ^( K Yr^-fj e (y>- x ) e (y< ) 
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(4.39) 

which reproduce the evolution equation for hi in master form 
dhi r 1 r l 

— = dywi(x\y)h!(y,T) - dyw 2 (y\x)h 1 (x,T). (4.40) 
dr Jo Jo 

The Kramers-Moyal expansion is derived as before, with some slight modifications. The 

result is obtained introducing an intermediate cutoff which is removed at the end. In this 

case we get 

dh x „ (YJ 2 3 6 n . (\ - x\\ . . s 

Notice that compared to the standard Fokker-Planck approximation, the boundary 
now generates a term on the left-hand-side of the equation proportional to q(x) which is 
absent in eq. (4.33). This and higher order approximations to the DGLAP equation can 
be studied systematically both analytically and numerically and it is possible to assess 
the validity of the approximation [25]. 



4.8 Helicities to LO 

As we have mentioned above, an interesting version of the usual DGLAP equation involves 
the helicity distributions. 

We start introducing [21] the DGLAP kernels for fixed helicities P ++ (z) = (P(z) + 
AP(z))/2 and P+-(z) = (P(z) - AP(z))/2 which will be used below. P(z) denotes 
(generically) the unpolarized kernels, while the AP(z) are the longitudinally polarized 
ones. These definitions, throughout the chapter, are meant to be expanded up to NLO, 
the order at which our numerical analysis holds. 

The equations, in the helicity basis, are 



dq+ (x) 



dt 



^(Pl q + (-)®q + (y) + Pl q 4-)®q-(y) 

2tt y y 

+p q + 9 + Q® 9+ (y) + P q + 9 -Q®9-(y)), 



dq_ 



x 



dt 



2ix 



(A 



,x 



,x. 



q+(y) + P++(-) ® q-(v) 

y y 
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4.8. Helicities to LO 
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Figure 4.2: LO kernels (qq and qg) in the helicity basis. 



+n g -(-) ® 9 + (v) + n 9 + (~) ® ff-(y)), 
= ® ? + (y) + ® My) 

at Zn y y 

+pr + (^) ® 9+ (y) + Pi 9 4^) ® g~m 

The nonsinglet (valence) analogue of this equation is also easy to write down 

dq + y(x) a s 



dt 

dq-y(x) 
dt 



-(P++(-) ® <z+,y(y) + P+_(-) ® g-,v(y)), 
2ti y y 

7^(P+-(j) ® ?+,v(y) + ® ?-,v(y))- 



(4.41) 



(4.42) 



where the g-ty = q± — q± are the valence components of fixed helicities. The kernels in 
this basis are given by 



3 (o) 

NS±,- 



P 



(0) 



99,H — 

3 (0) 
99>++ 
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n f x 2 
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99 
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Figure 4.3: NLO kernels (qq and qg) in the helicity basis. 
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1 



Pgq,— - C F 



N c [3x + - - 3-x' 

x 



x 



+ -- l- x-x 2 \ + p 5(l - x) 



(4.43) 



Taking linear combinations of these equations (adding and subtracting), one recovers 
the usual evolutions for unpolarized q{x) and longitudinally polarized Aq(x) distributions. 
We recall that the unpolarized distributions, the polarized and the transversely polarized 
qr{x) are related by 



q(x) = q + (x) + q_(x) = q +T (x) + q-r(x) 
Aq(x) = q+(x)-q4x) (4.44) 

at any Q of the evolution and, in particular, at the boundary of the evolution. 

Similar definition have been introduced for the gluon sector with G±(x) denoting 
the fixed helicities of the gluon distributions with Ag(x) = g+( x ) ~ 9-( x ) an( l di x ) = 
g + (x)+g-(x) being the corresponding longitudinal asymmetry and the unpolarized density 
respectively. 
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4.8. Helicities to LO 
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Figure 4.4: Nonsinglet kernels in the helicity basis. 
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Figure 4.5: LO kernels (gq and gg) in the helicity basis. 



Chapter 4. The kinetic interpretation of the DGLAP Equation 



105 



4.9 Summary of Positivity Arguments 

Let us recapitulate here the basic arguments [21] that are brought forward in order to 
prove the positivity of the evolution to NLO. 

If 

\AP(z)\ < P(z),z < 1 (4.45) 

then both kernels P ++ (z) and P + _(z) are positive as far as z < 1. 

The singular contributions at z — 1, which appear as subtraction terms in the evolution 
and which could, in principle, alter positivity, appear only in diagonal form, which means 
that they are only contained in P ++ , multiplied by the single functions q+(x) or q~(x) 

KM^rZtye^n... (4.40) 

*^ = £ m( £ )89+fe) + ... (4.47) 

^ = ^(f)®9 + (») + - ( 4 ' 48 > 

Let's focus just on the equation for q + (4.46). Rewriting the diagonal contribution as 
a master equation 

- = / dx' (w ++ (x\x')q + (x' , t) — w ++ (x'\x)q + (x, r)) dx' + ... (4.49) 
dt J 

in terms of a transition probability 

w ++ {x\y) = ^-P ++ (x/y)9(y > x) (4.50) 

which can be easily established to be positive, as we are going to show rigorously below, 
as far as all the remaining terms (the ellipses) are positive. We have performed a detailed 
numerical analysis to show the positivity of the contributions at x — 1. 

This last condition is also clearly satisfied, since the 5(1 — z) contributions appear only 
in P ++ and are diagonal in the helicity of the various flavours (q, g). For a rigorous proof 
of the positivity of the solutions of master equations we proceed as follows. 

Let q(x, r) be a positive distribution for r < r c and let us assume that it vanishes at 
r = t c , after which it turns negative. We also assume that the evolution of q(x, r) is of 
the form (4.7) with positive transition probabilities w(x\y) and w(y\x). Notice that since 
the function is continuous together with its first derivative and decreasing, continuity of 
its first derivative will require dg ^' T ' ) < at r = t c and in its neighbor. However, eq. (4.7) 
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Figure 4.6: NLO kernels (gq and gg) in the helicity basis. 



requires that at r = r c 



dq(x) 
dt 



= J dx'w(x\x')q(x',r) 



(4.51) 



which is positive, and we have a contradiction. We can picture the evolution in r of 
these functions as a family of curves getting support to smaller and smaller x- values as r 
grows and being almost vanishing at intermediate and large x values. We should mention 
that this proof does not require a complete probabilistic picture of the evolution, but just 
the positivity of the bulk part of the kernels, the positivity of the edge point subtractions 
and their diagonality in flavour. From Figs. (4.2) and (4.5) it is also evident that the 
leading order kernels are positive, together with the qg and qq (Fig. (4.3)) sectors. 

The edge point contributions, generating the "subtraction terms" in the master equa- 
tions for the "++" components of the kernels are positive, as is illustrated in 3 Tables 
included in Appendix A. There we have organized these terms in the form ~ CS(1 — x) 
with 

C = - log(l - A) A + B (4.52) 

with A and B being numerical coefficients depending on the number of flavours included 
in the kernels. Notice that the subtraction terms are always of the form (4.52), with 
the (diverging) logarithmic contribution (~ J A dzj{\ — z)) regulated by a cutoff. This 
divergence in the convolution cancels when these terms are combined with the divergence 
at x — 1 of the first term of the master equation for all the relevant components con- 
taining "+" distributions. It is crucial, however, to establish positivity of the evolution 
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Figure 4.7: The crossing of NLO kernels (gq and gg) in the helicity basis. 

of the helicities that the boundary conditions on the evolution \ Aq(x,Ql)\ < q(x,Ql) be 
satisfied. Initial conditions have this special property, in most of models, and the proof 
of positivity of all the distributions therefore holds at any Q. 

As we move to NLO, the pattern gets more complicated. In fact, from a numerical 
check, one can see that some NLO kernels turn to be negative, including the unpolarized 
kernels and the helicity kernels, while others (Fig. (4.4)) are positive. One can also notice 
the presence of a crossing of several helicity components in the gq and gg sectors (see Fig. 
(4.7)) at larger x values, while in the small-x region some components turn negative (Fig. 
(4.6)). There is no compelling proof of positivity, in this case, either than that coming 
from a direct numerical analysis. 

We have seen that master forms of evolution equations, for evolutions of all kinds, 
when found, can be used to establish positivity of the evolution itself. 

The requirements have been spelled out above and can be summarized in the following 
points: 1) diagonality of the decreasing terms, 2) initial positivity of the distributions, 3) 
positivity of the remaining (non diagonal) kernels. As we have also seen, some of these 
conditions are not satisfied by the NLO evolution. 

4.10 Results 

In regard to the input distributions we used, we refer the reader to Sec. 1.8. 

We show in Fig. 4.8 results for the evolution of the u + distribution at the initial scale 
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4.10. Results 



N f 
1 y j 


A 


B 


Q 
O 


1 9 ^09 


19 1 7^Q 


4 




1 D RQ94 


5 


9.3836 


9.2109 


6 


7.8103 


7.7249 



Table 4.1: Coefficients A and B for P, 



(0.632 GeV) and at two final scales, 100 GeV and 200 GeV respectively. The peaks at the 
various scales get lowered and become more pronounced toward the smaller x region as 
Q increases. In Fig. 4.9 u~ show an apparent steeper growth at small- x compared to u + . 
For the d distributions the situation is inverted, with d~ growing steeper compared to d + 
(Figs. 4.11 and 4.12 respectively). This apparent behaviour is resolved in Figs. (4.10) 
and (4.13) from which it is evident that both plus and minus components converge, at 
very small-x values, toward the same limit. 

The components s,c,t and b (Figs. 4.14-4.21) have been generated radiatively from 
vanishing initial conditions for final evolution scales of 100 and 200 GeV (s, c, b) and 200 
GeV (t). 

Both positive and negative components grow steadily at small-x and are negligible at 
larger x values. The distribution for the top quark (t) has been included for completeness. 
Given the smaller evolution interval the helicity distributions for heavier generations are 
suppressed compared to those of lighter flavours. Gluon helicities (Figs. 4.22, 4.23) are 
also enhanced at small-x, and show a similar growth. The fine difference between the 
quark u and d distributions are shown in Figs. 4.10 and 4.13. 

Finally in Figs. 4.24, 4.25 and 4.26 we plot simultaneously longitudinally polarized, 
unpolarized and helicity distributions for up quarks, down quarks and gluons at an in- 
termediate factorization scale of 100 GeV, relevant for experiments at RHIC. Notice that 
while u + and u~ are positive and their difference (Aw) is also positive, for down quarks 
the two helicity components are positive while their difference (Ad) is negative. Gluons, 
in the model studied here, have a positive longitudinal polarization, and their helicity 
components are also positive. The positive and negative gluon helicities are plotted in 
two separate figures, Figs. 4.22 and 4.23, while their difference, Ag(x) is shown in Fig. 
4.27. One can observe, at least in this model, a crossing at small-x in this distribution. 

We conclude that, at least for this set of boundary conditions, positivity of all the 
components holds to NLO, as expected. 

We also report three tables (4.1 - 4.3) illustrating the (positive) numerical values of 
the contributions coming from the subtraction terms in the NLO kernels. Coefficients A 
and B refer to the subtraction terms — log(l — A)A + B as explained in the section above. 
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Table 4.2 
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48.4555 
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24.0579 
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20.7245 
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40.4555 


17.3912 



qq,++ 



Table 4.3: Coefficients A and B as in Table 4.1 for Pgl] ++ 

4.11 Regularizations 

The "+" plus form of the kernels and all the other forms introduced before, obtained by 
separating the contributions from the edge-point (x — 1) from those coming from the 
bulk (0 < x < 1) are all equivalent, as we are going to show, with the understanding 
that a linear (unique) cutoff is used to regulate the divergences both at x=0 and at x=l. 
We focus here on the two possible sources of singularity, i.e. on P qq and on the P gg 
contributions, which require some attention. Let's start from the P qq case. We recall that 
"+" plus distributions are defined as 



1 



6(1 — x — A) 



5(1 



(1 — x)+ 1 — x 

with A being a cutoff for the edge-point contribution. 
We will be using the relations 



x) 



l-A 

o T 



dz 



(4.53) 



Jx y 



1 dz 
o 1 — z 



1 dz 



dy 



f(y)g(x/y) 



x 1 — Z 



log(l - x) 



ix y Jx y 

Using the expressions above it is easy to obtain 



f(x/y)g(y)- 



(4.54) 



(1 - x)+ Jx y (1 -x/y)+ 

Jx y 1 — x/y^^ Jx y ^ ^ Jo z 
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1 dy yfjy) - xf(x) 

= / + log(l - x)f(x) 

Jx y y-x 

(4.55) 

which is eq. (4.18). If we remove the "+" distributions and adopt (implicitly) a cutoff 
regularization, we need special care. In the probabilistic version of the kernel, the handling 
of P qq <S> q is rather straightforward 



l-x/y) q{y) 



x/y 

- c r£y-^o<Th- 1 -*) 

(4.56) 

and using eqs. (4.54) we easily obtain 

Jx y y — x 

-C F + x/y) q(y) + ^-C F q{x). (4.57) 

Jx y i 

Now consider the convolution P gg <g> g(x) in the Collins-Qiu form. We get 

p M 0, w = 2C ,/;| T -L 7 - 9(!/ ) 

+2Ca /; * ( lM1 _ I/v) + _l _ 2 ) 9(B) _ sW /; ± a (i + 

-«?(x) jT 1 2C A (z(l -z)-2)- n f g(x) £ dz 1 - (z 2 + (1 - zf) . (4.58) 



1 

T 



There are some terms in the expression above that require some care. The appropriate 
regularization is 

r 1 dz f 1 dz T/A . r 1 dz r 1 ~ A dz 

/ — + / i ► / A = / — + / . 4.59 

JO 2 JO 1 — Z J A 2 JO 1 — 2 



Observe also that 

r 1 dz r 1 - A dz 

— = z = -logA. 4.60 

JA Z JO 1 — Z 

Notice that in this regularization the singularity of 1/z at z = is traded for a singularity 
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at z=l in 1/(1 — z). It is then rather straightforward to show that 

dyyg(y) - xg(x) 



P gg <g> g(x) = 2C A f 1 



y y-x 

dy 

y 



+ 2C A \og(l-x)g(x) 



1 



-2C A /'* — {, ;! I -,•/,/) — -2 



A) 



<?(*). (4.61) 



A final comment is due for the form of the sum rule 



Si dz(z - \) l99 = 



(4.62) 



that we need to check with the regularization given above. 

The strategy to handle this expression is the same as before. We extract all the 1/z 
and 1/(1 — z) integration terms and use 



r 1 dz r 1 dz r 1 
Jo z Jo 1 — z J A 



1 dz fi-A dz 



1 — Z J A Z JO 1 — Z 

= -logA + logA = 



(4.63) 



to eliminate the singularities at the boundaries x — 0, 1 and verify eq. (4.62). 
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Figure 4.8: Evolution of u + versus x at various Q values. 
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Figure 4.11: Evolution of d + versus x at various Q values. 



x 



0.35 
0.3 

0.25 
0.2 

0.15 
0.1 

0.05 




1 1 1 


Q=0.632GeV 




Q=100GeV 




Q=200GeV - 


/ "'-\ 
/ "***- * 

/ 

-/ ***'.* 




1 1 1 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x 



Figure 4.12: Evolution of d versus x at various Q values. 
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Figure 4.13: Small-x behaviour of at 100 GeV. 
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Figure 4.14: Evolution of s + versus x at various Q values. 



Chapter 4. The kinetic interpretation of the DGLAP Equation 



115 



0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 




Q=0.632 GeV 
Q=100 GeV 
Q=200 GeV 



0.14 
0.12 
0.1 
0.08 



£ 0.06 



0.04 



0.02 



0=0.632 GeV 
Q=100 GeV 
Q=200 GeV 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x 

Figure 4.15: Evolution of s~ versus x at various Q values. 
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Figure 4.16: Evolution of c + versus x at various Q values. 
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Figure 4.17: Evolution of c~ versus x at various Q values. 
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Figure 4.18: Evolution of b + versus x at various Q values. 
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Figure 4.19: Evolution of b versus x at various Q values. 



4.12 Conclusions 

We have discussed in detail some of the main features of the probabilistic approach to 
the DGLAP evolution in the helicity basis. Numerical results for the evolution of all the 
helicities have been provided, using a special algorithm, based in x-space. We have also 
illustrated some of the essential differences between the standard distributional form of 
the kernels and their probabilistic version, clarifying some issues connected to their regu- 
larization. Then we have turned to the probabilistic picture, stressing on the connection 
between the random walk approach to parton diffusion in x-space and the master form 
of the DGLAP equation. The link between the two descriptions has been discussed es- 
pecially in the context of the Kramers-Moyal expansion. A Fokker-Planck approximation 
to the expansion has also been presented which may turn useful for the study of formal 
properties of the probabilistic evolution. We have also seen that positivity of the helicity 
distributions, to NLO, requires a numerical analysis, as already hinted in [21]. Our study 
also validates the use of a very fast evolution algorithm, alternative to other standard 
algorithms based on Mellin algorithms, whose advantage is especially in the analysis of 
the evolution of nonforward parton distributions, as we will show elsewhere. 
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Figure 4.20: Evolution of t + versus x at various Q values. 
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Figure 4.21: Evolution oft versus x at various Q values. 
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Figure 4.22: Evolution of g + versus x at various Q values. 
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Figure 4.23: Evolution of g versus x at various Q values. 
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Figure 4.26: Various kinds of gluon distributions at Q = 100 GeV. 
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Figure 4.27: Small-x behaviour of Ag at 100 GeV. 
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Chapter 5 



An x-space analysis of evolution 
equations: Soffer's inequality and the 
nonforward evolution 

We analyze the use of algorithms based in x-space for the solution of renormalization group 
equations of DGLAP-type and test their consistency by studying bounds among partons 
distributions - in our specific case Soffer's inequality and the perturbative behaviour of 
the nucleon tensor charge - to next-to-leading order in QCD. 

We also comment on the (kinetic) proof of positivity of the evolution of h\ , using a 
kinetic analogy, along the lines of the previous chapter, and illustrate the extension of the 
algorithm to the evolution of generalized parton distributions. We prove positivity of the 
nonforward evolution in a special case and illustrate a Fokker-Planck approximation to 
it. 

5.1 Introduction 

One of the most fascinating aspects of the structure of the nucleon is the study of the 
distribution of spin among its constituents, a topic of remarkable conceptual complexity 
which has gained a lot of attention in recent years. This study is entirely based on 
the classification and on the phenomenological modeling of all the leading-twist parton 
distributions, used as building blocks for further investigations in hadronic physics. 

There are various theoretical ways to gather information on these non-local matrix 
elements. One among the various possibilities is to discover sum rules connecting mo- 
ments of these distributions to other fundamental observables. Another possibility is to 
discover bounds - or inequalities - among them and use these results in the process of 
their modeling. There are various bounds that can be studied, particularly in the context 
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of the new generalized parton dynamics typical of the skewed distributions [82, 106]. All 
these relations can be analyzed in perturbation theory and studied using the Renormal- 
ization Group (RG), although a complete description of their perturbative dynamics is 
still missing. This study, we believe, may require considerable theoretical effort since it 
involves a global understanding both of the (older) forward (DGLAP) dynamics and of 
the generalized new dynamics encoded in the skewed distributions. 

In this context, a program aimed at the study of various bounds in perturbation the- 
ory using primarily a parton dynamics in x-space has been outlined [25]. This requires 
accurate algorithms to solve the equations up to next-to-leading order (NLO). Also, un- 
derlying this type of description is, in many cases, a probabilistic approach [22] which 
has some interesting consequences worth of a closer look . In fact, the DGLAP equation, 
viewed as a probabilistic process, can be rewritten in a master form which is at the root of 
some interesting formal developments. In particular, a wide set of results, available from 
the theory of stochastic processes, find their way in the study of the evolution. We have 
elaborated on this issue in previous work [25] and proposed a Kramers-Moyal expansion of 
the DGLAP equation as an alternative way to describe its dynamics. Here, this analysis 
will be extended to the case of the nonforward evolution. 

With these objectives in mind, in this study we test x-space algorithms up to NLO and 
verify their accuracy using a stringent test: Soffer's inequality. As usual, we are bound to 
work with specific models of initial conditions. The implementations on which our analysis 
are based are general, with a varying flavour number n/ at any threshold of intermediate 
quark mass in the evolution. Here, we address Soffer's inequality using an approach based 
on the notion of "superdistributions" [21], which are constructs designed to have a simple 
(positive) evolution thanks to the existence of an underlying master form [22, 25]. The 
original motivation for using such a master form (also termed kinetic or probabilistic) to 
prove positivity has been presented in [21], while further extensions of these arguments 
have been presented in [25]. In a final section we propose the extension of the evolution 
algorithm to the case of the skewed distributions, and illustrate its implementation in the 
nonsinglet case. As for the forward case, numerical tests of the inequality are performed 
for two different models. We show that even starting from a saturated inequality at 
the lowest evolution scale, the various models differ significantly even for a moderate 
final factorization scale of Q = 100 GeV. Finally, we illustrate in another application the 
evolution of the tensor charge and show that, in the models considered, differences in the 
prediction of the tensor charge are large. 
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5.2 Prelude to x-space: A Simple Proof of Positivity of 
hi to NLO 

There are some nice features of the parton dynamics, at least in the leading logarithmic 
approximation (LO), when viewed in x-space, once a suitable "master form" of the parton 
evolution equations is identified. 

The existence of such a master form, as firstly shown by Teryaev, is a special feature 
of the evolution equation itself. The topic has been addressed before in LO [21] and 
reanalyzed in more detail in [25] where, starting from a kinetic interpretation of the evo- 
lution, a differential equation obtained from the Kramers-Moyal expansion of the DGLAP 
equation has also been proposed. 

The arguments of refs. [21, 25] are built around a form of the evolution equation 
which has a simple kinetic interpretation and is written in terms of transition probabilities 
constructed from the kernels. 

The strategy used, at least in leading order, to demonstrate the positivity of the LO 
evolution for special combinations of parton distributions Q± [21], to be defined below, 
or the NLO evolution for hi, which we are going to address, is based on some results of 
ref.[21], briefly reviewed here, in order to be self-contained. 

A master equation is typically given by 

d f 

f(x,r)= / dx' (w(x\x')f(x', t) — w(x'\x)f(x, r)) dx' (5.1) 



dr 

and if through some manipulations, a DGLAP equation 
with kernels P(x), is rewritten in such a way to resemble eq. (5.1) 

U{ X , r) = f 1 dyP (A^il_r^yp (i) *m (5 . 3) 

dr Jx \U J U Jo y \xj x 

with a (positive) transition probability 

w(x\y) = ^P(x/y) 9 -^l (5.4) 
2tt y 

then positivity of the evolution is established. 

For equations of nonsinglet type, such as those evolving = q — q, the valence 
quark distribution, or hi, the transverse spin distribution, this rewriting of the equation 
is possible, at least in LO. NLO proofs are, in general, impossible to construct by this 
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method, since the kernels turn out, in many cases, to be negative. The only possible 
proof, in these cases, is just a numerical one, for suitable (positive) boundary conditions 
observed by the initial form of the parton distributions. Positivity of the evolution is then 
a result of a non obvious interplay between the various contributions to the kernels in 
various regions in x-space. 

In order to discuss the probabilistic version of the DGLAP equation it is convenient to 
separate the bulk contributions of the kernels (x < 1) from the edge point contributions 
at x = 1. For this purpose we recall that the structure of the kernels is, in general, given 
by 

P(z) = P(z) - 5(1 - z) [ P(z)dz, (5.5) 

Jo 

where the bulk contributions (z < 1) and the edge point contributions (~ S(z — 1)) have 
been explicitly separated. We focus on the transverse spin distributions as an example. 
With these prerequisites, proving the LO and NLO positivity of the transverse spin dis- 
tributions is quite straightforward, but requires a numerical inspection of the transverse 
kernels. Since the evolutions for Axq^ = h\ are purely nonsinglet, diagonality in flavour 
of the subtraction terms (~ Jq w(y\x)q(x, r)) is satisfied, while the edge-point subtrac- 
tions can be tested to be positive numerically. We illustrate the explicit construction of 
the master equation for hi in LO, since extensions to NLO of this construction are rather 
straightforward. 

In this case the LO kernel is given by 



A T P£\x) = C F 



(T^xT + - 2 + >-^ 



(5.6) 



and by some simple manipulations we can rewrite the corresponding evolution equation 
in a suitable master form. That this is possible is an elementary fact since the subtraction 
terms can be written as integrals of a positive function. For instance, a possibility is to 
choose the transition probabilities 



wi[x\y] = — 



- 2 % > x)0(y < 1) 



w 2 [y\x] = 



X 



l-y/x 2 / 



9{y > -x)9{y < 0) 



(5.7) 



which reproduce the evolution equation for h\ in master form 



dhi 



= [ dywi{x\y)hi{y,r) - f dyw 2 (y\x)h 1 (x,r). 
Jo Jo 



(5.8) 
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fir 


A 


B 


Q 
O 


1 9 ^09 


19 1 7^Q 


4 




1 D RQ94 


5 


9.3836 


9.2109 


6 


7.8103 


7.7249 



Table 5.1: Coefficients A and B for various flavor numbers, to NLO for AxPqq^ 

The NLO proof of positivity is also rather straightforward. For this purpose we have 
analyzed numerically the behaviour of the NLO kernels both in their bulk region and at 
the edge-point. We show in Table 1 of Appendix B results for the edge point contributions 
to NLO for both of the A T P±^ components, which are numerically the same. There we 
have organized these terms in the form ~ C5(l — x) with 

C = -log(l - A)A + B, (5.9) 

with A and B being numerical coefficients depending on the number of flavours included 
in the kernels. The (diverging) logarithmic contribution (~ J A dz/(l — z)) have been 
regulated by a cutoff. This divergence in the convolution cancels when these terms are 
combined with the divergence at x = 1 of the first term of the master equation (5.8) for 
all the relevant components containing "+" distributions. As for the bulk contributions 
(x < 1), positivity up to NLO of the transverse kernels is shown numerically in Fig. 
(5.1). All the conditions of positivity are therefore satisfied and therefore the A T ±q 
distributions evolve positively up to NLO. The existence of a master form of the equation 
is then guaranteed. 

Notice that the NLO positivity of Ax±q implies positivity of the nucleon tensor charge 

[81] 

8q = j*dx(h\{x)-hf x) ) (5.10) 

for each separate flavour for positive initial conditions. As we have just shown, this proof 
of positivity is very short, as far as one can check numerically that both components of 
eq.(5.8) are positive. 



5.3 Soffer's inequality 

Numerical tests of Soffer's inequality can be performed either in moment space or, as we 
are going to illustrate in the next section, directly in x-space, using suitable algorithms 
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to capture the perturbative nature of the evolution. We recall that Soffer's inequality 

Ihix)] <q + {x) (5.11) 

sets a bound on the transverse spin distribution hi(x) in terms of the components of 
the positive helicity component of the quarks, for a given flavour. An original proof of 
Soffer's inequality in LO has been discussed in ref.[14], while in [21] an alternative proof 
was presented, based on a kinetic interpretation of the evolution equations. 

We recall that hi, also denoted by the symbol 

A T q(x, Q 2 ) = ql(x, Q 2 ) - J(x, Q 2 ), (5.12) 

has the property of being purely nonsinglet and of appearing at leading twist. It is 
identifiable in transversely polarized hadron-hadron collisions and not in Deep Inelastic 
Scattering (from now on we will omit sometimes the x-dependence in the kernels and 
in the distributions when obvious). In the following we will use interchangeably the 
notations hi = h\ and At<? to denote the transverse asymmetries. We introduce also the 
combinations 

A T (q + q) = hl + hl 
A T q^=A T (q-q) = h\-h\ 

A T gW = £)A r ( ft + fc) 

i 

(5.13) 

where we sum over the flavor index (i), and we have introduced singlet and nonsinglet 
contributions for distributions of fixed helicities 

q+ ] = J2(<i+i + <i+i) 

i 

Q+ ] = Q+i - g+i = S. 

(5.14) 

In our analysis we solve all the equations in the helicity basis and reconstruct the various 
helicities after separating singlet and nonsinglet sectors. We mention that the nonsinglet 
sector is now given by a set of 2 equations, each involving ± helicities and the singlet 
sector is given by a 4-by-4 matrix. 

In the singlet sector we have 
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dG + (x) _ a s Gq ( +) Gq ( +) 



dt 2tt 



+P GG ®G' + + P GG ®GL), 

+P GG ®G + + P GG ®GL). (5.15) 
while the nonsinglet (valence) analogue of this equation is also easy to write down 

^r 1 = £ (p +- q + ] + p ++ ^ )} - (5 - 16) 

Above, i is the flavor index, (±) indicate q ± g components and the lower subscript ± 
stands for the helicity. 

Similarly to the unpolarized case the flavour reconstruction is done by adding two 
additional equations for each flavour in the helicity ± 

X±,i = ~ -q ( ± ] (5-17) 



whose evolution is given by 



dt 27T 1 ++ 



dx+ ~ i{x) - as (p:?®x +i +p??®x^, 



(5.18) 



The reconstruction of the various contributions in flavour space for the two helicities 
is finally done using the linear combinations 

q±i = \ (q±i + x±* + ^<?i +) ) • (5.19) 
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We will be needing these equations below when we present a proof of positivity up 
to LO, and we will thereafter proceed with a NLO implementation of these and other 
evolution equations. For this we will be needing some more notations. 

We recall that the following relations are also true to all orders 

P(x) = ±(P ++ (x) + P + -(x)) 
= l(P-(x) + P_ + (x)) 

between polarized and unpolarized (P) kernels and 

P ++ (x)=P^(x), P_ + (x) = P + _(x) (5.20) 

relating unpolarized kernels to longitudinally polarized ones. Generically, the kernels of 
various type are expanded up to NLO as 

P(x) = ^\x) + ^) 2 P^(X), (5.21) 
and specifically, in the transverse case we have 



A T P« ± = A T P«±A T P«, (5.22) 

(5.23) 



with the corresponding evolution equations 

d 



& T q±(Q 2 ) = A T P qq , ± (a s (Q 2 )) ® A T q ± (Q 2 ) . (5.24) 



cilnQ 2 

We also recall that the kernels in the helicity basis in LO are given by 



p(0) 

^NS±,++ 


p(0) 


- p(0) 

qq 


p(0) 


p(0) 


= 


p(0) 


= n f x 2 




Pqg,+- 


= Pqg,-+ 


= n f (x 

X 


Pgq,++ 


— P 

r gq, — 


p(0) 


. p(0) 


= N C ( 



X J _i_ X 



1-x-x 2 ) +/3o6(l-x) 



P^ + _ = N c ( 3x + \ - 3 - x 2 ) . (5.25) 



x 
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An inequality, such as Soffer's inequality, can be stated as positivity condition for 
suitable linear combinations of parton distributions [21] and this condition can be analyzed 
- as we have just shown for the hi case - in a most direct way using the master form. 

For this purpose consider the linear valence combinations 

Q+ = <?+ + h 
Q = q+- hi 

(5.26) 

which are termed "superdistributions" in ref.[21]. Notice that a proof of positivity of 
the Q distributions is equivalent to verify Soffer's inequality. However, given the mixing 
of singlet and nonsinglet sectors, the analysis of the master form is, in this case, more 
complex. As we have just mentioned, what can spoil the proof of positivity, in general, is 
the negativity of the kernels to higher order. We anticipate here the result that we will 
illustrate below where we show that a LO proof of the positivity of the evolution for Q 
can be established using kinetic arguments, being the kernels are positive at this order. 
However we find that the NLO kernels do not satisfy this condition. In any case, let's see 
how the identification of such master form proceeds in general. We find useful to illustrate 
the result using the separation between singlet and nonsinglet sectors. In this case we 
introduce the combinations 



Q<-> = ? H±M-» 



with h ( f s A T 9«. 

Differentiating these two linear combinations (5.27) we get 



d\og(Q 2 ) 



which can be rewritten as 



(5.27) 



(5.28) 



- 5 (Pil + #>) OP + i (4;> - pV) qL-> + pi-'-'-' 



q- 



rflog(Q 2 ) 2 V ++ J 1 2 



rflog(Q 2 ) 
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Figure 5.1: Plot of the transverse kernels. 



with ) = P NS being the nonsinglet (NS) kernel. 
At this point we define the linear combinations 



P?± = 2 ( P ++ ± p r) 



and rewrite the equations above as 



P? + Q i+ + P?-Q i - + n q -Qi- 



d\og(Q 2 ) 

rfio g (g 2 ) 



(5.29) 



(5.30) 



(5.31) 



where we have reintroduced i as a flavour index. From this form of the equations it is 
easy to establish the leading order positivity of the evolution, after checking the positivity 
of the kernel and the existence of a master form. 

The second nonsinglet sector is defined via the variables 



Xi± - Qi± Qi± 



(5.32) 
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which evolve as nonsinglets and the two additional distributions 

Q xli± = X ,+ ±/4 (+) . (5-33) 
Also in this case we introduce the kernels 

P?l = l(P ++ ±A T P^) (5.34) 

to obtain the evolutions 

dQxi+ _ pQxn _l r>Qxr\ _l mi 



dlog(Q 2 ) 



(5.35) 



For the singlet sector, we simply define Q+~^ = q^ + \ and the corresponding evolution 
is similar to the singlet equation of the helicity basis. Using the equations above, the 
distributions Qj± are then reconstructed as 




Q,± = o Off + Qx4 + —Q+ ( 5 - 36 ) 



and result positive for any flavour if the addends are positive as well. However, as we have 
just mentioned, positivity of all the kernels introduced above is easy to check numerically 
to LO, together with their diagonality in flavour which guarantees the existence of a 
master form. 

As an example, consider the LO evolution of Q±. The proof of positivity is a simple 
consequence of the structure of eq. (5.31). In fact the edge-point contributions appear 
only in -P++, i.e. they are diagonal in the evolution of Q±. The inhomogeneous terms on 
the right hand side of (5.31), proportional to q_ are are harmless, since the P + _ kernel has 
no edge-point contributions. Therefore under 1) diagonality in flavour of the subtraction 
terms and 2) positivity of first and second term (transition probabilities) we can have 
positivity of the evolution. A refined arguments to support this claim has been presented 
in [25]. 

This construction is not valid to NLO. In fact, while the features of flavour diagonality 
of the master equation are satisfied, the transition probabilities w(x,y) are not positive 
in the whole x, y range. The existence of a crossing from positive to negative values in 
P+ + can, in fact, be established quite easily using a numerical analysis. We illustrate in 
Figs. (5.2) and (5.3) plots of the Q kernels at LO and NLO, showing that, at NLO, the 
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Figure 5.2: Plot of the LO kernels for the Q distributions 



requirement of positivity of some components is violated. The limitations of this sort of 
proofs -based on kinetic arguments- are strictly linked to the positivity of the transition 
probabilities once a master form of the equation is identified. 



5.4 An x-space Expansion 

We have seen that NLO proofs of positivity, can be -at least partially- obtained only 
for suitable sets of boundary conditions. To this purpose, we choose to investigate the 
numerical behaviour of the solution using x-space based algorithms which need to be 
tested up to NLO. 

Our study validates a method which can be used to solve evolution equations with 
accuracy in leading and in next-to-leading order. The method is entirely based on an 
expansion [109] used in the context of spin physics [75] and in supersymmetry [39]. An 
interesting feature of the expansion, once combined with Soffer's inequality, is to gen- 
erate an infinite set of relations among the scale invariant coefficients (A n ,B n ) which 
characterize it. 

In this approach, the NLO expansion of the distributions in the DGLAP equation is 
generically given by 

f(x ' > = S — log [«Qf)) + a(Q > £ — log (5 - 37) 
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Figure 5.3: Plot of the NLO kernels for the Q distributions, showing a negative behaviour 
at large x 



where, to simplify the notation, we assume a short-hand matrix notation for all the 
convolution products. Therefore f(x,Q 2 ) stands for a vector having as components any 
of the helicities of the various flavours (Q±, q±, G±). 

If we introduce Rossi's expansion for hi, q+, and the linear combinations Q± (in short 
form) 



q ± ~ (A?,B?) 

(5.38) 

we easily get the inequalities 

(-1)"(^++^) >0 (5.39) 

and 

(_!)« (Af - A^j > (5.40) 

valid to leading order, which we can check numerically. Notice that the signature factor 
has to be included due to the alternation in sign of the expansion. To next to leading 
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order we obtain 

(-ir +1 (Al+(x) + a(Q*)B«+(x)) < (-l) n (A h n (x) + a(Q*)B h n (x)) 

< (-l) n (A^(x) + a(Q 2 )B^(x)) (5.41) 

valid for n > 1, obtained after identification of the corresponding logarithmic powers 
\og(a(Q 2 )) at any Q. In general, one can assume a saturation of the inequality at the 
initial evolution scale 

Q_(x, Ql) = h^x, Ql) - ^q + (x, Ql) = 0. (5.42) 

This initial condition has been evolved in Q solving the equations for the Q± distributions 
to NLO. 

5.5 Nonforward Extensions 

In this section we finally discuss the nonforward extension of the evolution algorithm. In 
the case of nonforward distributions a second scaling parameter ( controls the asymmetry 
between the initial and the final nucleon momentum in the deeply virtual limit of nucleon 
Compton scattering. The solution of the evolution equations, in this case, are known in 
operatorial form. Single and double parton distributions are obtained sandwiching the 
operatorial solution with 4 possible types of initial/final states < >, < p|--|0 > 

,< p'\---\p >, corresponding, respectively, to the case of diagonal parton distributions, 
distribution amplitudes and, in the latter case, skewed and double parton distributions 
[106]. Here we will simply analyze the nonsinglet case and discuss the extension of the 
forward algorithm to this more general case. Therefore, given the off-forward distributions 
H q (x,£), in Ji's notation, we set up the expansion 

= £ H log Uwi) j + a(Q » £ — iti— log Uwi) j ■ (6 - 43) 

which is the natural extension of the forward algorithm discussed in the previous sections. 
We recall that in the light-cone gauge H(x,£) is defined as 

# ? (*,£,A 2 )) = i/^e-^-(P'|^(0,^,0 ± )i 7 + ^(0,^,0 ± )|P) (5.44) 

with A = P' - P, P+ = 1/2 (P + P) [82] (symmetric choice) and £P = 1/2 A+. 

This distribution describes for x > £ and x < — £ the DGLAP-type region for the 
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quark and the antiquark distribution respectively, and the ERBL [52] (see also [44] for 
an overview) distribution amplitude for — £ < x < £. In the following we will omit the A 
dependence from H q . 

Again, once we insert the ansatz (5.43) into the evolution equations we obtain an 
infinite set of recursion relations which we can solve numerically. In LO, it is rather simple 
to relate the Gegenbauer moments of the skewed distributions and those of the generalized 
scaling coefficients A n . We recall that in the nonforward evolution, the multiplicatively 
renormalizable operators appearing in the light cone expansion are given in terms of 
Gegenbauer polynomials [106]. The Gegenbauer moments of the coefficients A n of our 
expansion (5.43) can be easily related to those of the off-forward distribution 

Cn(£,Q 2 ) = C f_Cl /2 {z/i)H{z^Q 2 )dz. (5.45) 
The evolution of these moments is rather simple 

c»(C,Q 2 ) = c»(C,0g)[^J (5-46) 

with 

1 2±J 1 




2 (n + l)(n + 2) ' %] j {bA1) 

being the nonsinglet anomalous dimensions. If we define the Gegenbauer moments of our 
expansion 

4 n) (£,Q 2 ) =c£ i C$ 2 (z/Z)H(zZ,Q 2 )dz (5.48) 
we can relate the moments of the two expansions as 

4°(0 = C„(C,Q^g)*- (5.49) 

Notice that expansions similar to (5.43) hold also for other choices of kinematical variables, 
such as those defining the nonforward distributions [106], where the t-channel longitudinal 
momentum exchange A + is related to the longitudinal momentum of the incoming nucleon 
as A = £P. We recall that H q (x.£) as defined in [82] can be mapped into two independent 
distributions F q (X, () and T q {X, () through the mappings [72] 



X 1 = 

x 2 = 



(*i + Q 
(i + O 

£ -3?2 
(1+0 

C/(2-C) 
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1 



l-C/2"' ! 

-1 



C/2 



(5.50) 



in which the interval — 1 < x < 1 is split into two coverings, partially overlapping (for 
— i < x < £, or ERBL region) in terms of the two variables — £ < Xi < 1 (0 < X x < 1) 
and —1 < X2 < £ (0 < X2 < 1). In this new parameterization, the momentum fraction 
carried by the emitted quark is X, as in the case of ordinary distributions, where it is 
parametrized by Bjorken x. For definitess, we focus here on the DGLAP-like (X > () 
region of the nonsinglet evolution. The nonsinglet kernel is given in this case by (x = X) 



7r \y — x 



1 + 



XX 



yy' 



f 1 1 ■ 

5(x-y) / dz— 
Jo 1 



(5.51) 



we introduce a LO ansatz 



(5.52) 



k=0 



and insert it into the evolution of this region to obtain the very simple recursion relations 

An +1 (X,0 = 



_1 Cf f 1 dy yMy,0-xMX,0 _ ^ Cp f 1 dy(x-Q (yA n (x,() -Ufe,Q) 

Po F Jx y y-X p F Jx y(y- () y-X 



--C F A n (X,() 

Po 



3 (i-xni-x/o 

2 l-c 



(5.53) 



The recursion relations can be easily reduced to a weighted sum of contributions in which 
£ is a spectator parameter. Here we will not make a complete implementation, but we 
will illustrate in an appendix the general strategy to be followed. There we show a very 
accurate analytical method to evaluate the logarithms generated by the expansion without 
having to rely on brute-force computations. 



5.6 Positivity of the nonsinglet Evolution 

Positivity of the nonsinglet evolution is a simple consequence of the master- form associated 
to the nonforward kernel (5.51). As we have already emphasized above, positivity of the 
initial conditions are sufficient to guarantee a positivity of the solution at any scale Q. 
The master-form of the equation allows to reinterpret the parton dynamics as a random 
walk biased toward small-x values as r = \og(Q 2 ) increases. 
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In the nonforward case the identification of a transition probability for the random 
walk [25] associated with the evolution of the parton distribution is obtained via the 
nonforward transition probability 



w c (x\y) 
w' c (y\x) 



a 1 
— u F 

7r y — x 

-c F t + y \ o(y < x) 



i x (x-C) 

y y-C 



9(y > x) 



n x 2 (x — y) 
and the corresponding master equation is given by 



(5.54) 



dT q 
~dT 



rl rx 

/ dyw c (x|y)^" ? (y,C,r) - / dy w'Ay\x)F q (x, (, r), (5.55) 

Jx JO 



that can be re-expressed in a form which is a simple generalization of the formula for the 
forward evolution [25]. 



"1 rx 

S i ( )2 ^ dyw c (x\y)F q (y, (, r) - ^ dy w'^y^F^x, (, r) 

ra(x) r—x 

= -J dyw ( (x + y\x)*J r q (x,(,T) + J^ dyw[(x + y\x)T q {x, C, ^.56) 



where a Moyal-like product appears 



W((x + y\x) * J- q {x, (, t) = w^(x + y\x)e 



(5.57) 



and a(x) — x — 1. A Kramers-Moyal expansion of the equation allows to generate a 
differential equation of infinite order with a parametric dependence on ( 



dT q 



rO r—x 

= / dyw c (x + y\x)J 7 q (x,(,r)+ dyw'^x + y^F^xX^) 

GElOgC^/ J a(x) JO 

~ E / dy K -^-d x n ( W< (x + y\x)F q (x, (, r)) . (5.58) 
n=1 J0 n. 



We define 



rO r—x 

d (x,() = / dyw c (x + y\x)J 7 q (x,(,T)+ dy w' ( (x + y^F^x, C, t) 

Ja(a;) JO 

a n (x,C) = / dyy n w c (x + y|x)^(x,C,r) 

« 

d n (x,C) = r (X) dyy n d x n (w c (x + y\x)F q (x,(,r)) n = l,2,... (5.59) 
Jo 
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If we arrest the expansion at the first two terms (n = 1, 2) we are able to derive an 
approximate equation describing the dynamics of partons for non-diagonal transitions. 
The procedure is a slight generalization of the method presented in [25], to which we refer 
for further details. For this purpose we use the identities 



ai(x,Q = d x ai(x, — a(x)d x a(x)w^(x + a(a;)|a;)jF (? (x, (, r) 
a 2 (x,() = dla 2 (x, () - 2a(x)(d x a(x)) 2 w <: (x + a(x)\x)J 7 q (x, (, r) 
-a(x) 2 d x a(x)d x (w c (x + a(x)\x)F q (x, (, r)) 
-a 2 (x)d x a(x)d x (w c (x + y\x)F q (x, (, r)) \ y=a (x)- (5.60) 

which allow to compute the first few coefficients of the expansion. Using these relations, 
the Fokker-Planck approximation to this equation can be worked out explicitly. We omit 
details on the derivation which is non obvious since particular care is needed to regulate 
the (canceling) divergences and just quote the result. 
A lengthy computation gives 



dr 



a n ( x °- 3 ^o,-i , W/ * \ 

-, Cf {{x-cr + + Xo '° ) ^ c ' r) 

+ fa ((^ + V^o) d ^ Cr) + ^(f=^^ (x ' c ' r) 



(5.61) 



where we have defined 

-1 + x) 3 (17x 3 - C 2 (3 + 4C) + 3x( (3 + 5C) - 3x 2 (3 + 7C))) 



^0-3 
#0,-1 



12 x 3 

-29x 4 - 3 + x 2 (-1 + C) + 2C - 2x (1 + 30 + x 3 (12 + 230 



3a; 3 



1 3 rtl (1-x) 
^o,o = 4+ — -- + 21og v 



Xi-i 



2x 2 x x 
- ((-1 + Qx - 15x 2 + 14x 3 ) (x - 0) 
3x~ 2 



1 5x , 23x 4 7( 3C 5x( 
Xl >- 3 = 2-Y + 5X - — + T-4x- + — 

2> 131a: 3 C 5C 2 ( 2 ( 2 >2 39x 2 ( 2 - C 3 8a< 3 



-((-l + x y( x -(y ( 3 + 23:r 2 + 4( - 2* (7 + 8(_ 
X2,-s - • (5 - 62) 



This equation and all the equations obtained by arresting the Kramers-Moyal expan- 
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sion to higher order provide a complementary description of the nonforward dynamics in 
the DGLAP region, at least in the nonsinglet case. Moving to higher order is straightfor- 
ward although the results are slightly lengthier. A full-fledged study of these equations is 
under way and we expect that the DGLAP dynamics is reobtained - directly from these 
equations - as the order of the approximation increases. 

5.7 Model comparisons, saturation and the tensor charge 

In this last section we discuss some implementations of our methods to the standard 
(forward) evolution by doing a NLO model comparisons both in the analysis of Soffer's 
inequality and for the evolution of the tensor charge. We have selected two models, 
motivated quite independently and we have compared the predicted evolution of the 
Soffer bound at an accessable final evolution scale around 100 GeV for the light quarks 
and around 200 GeV for the heavier generations. At this point we recall that in order 
to generate suitable initial conditions for the analysis of Soffer's inequality, one needs an 
ansatz in order to quantify the difference between its left-hand side and right-hand side 
at its initial value. 

The well known strategy to build reasonable initial conditions for the transverse spin 
distribution consists in generating polarized distributions (starting from the unpolarized 
ones) and then saturate the inequality at some lowest scale, which is the approach we 
have followed for all the models that we have implemented. 

The first model we have used is the GRV-GRSV, that we have described in Sec. 1.8. 

In the implementation of the second model (GGR model) we have used as input 
distributions in the unpolarized case the CTEQ4 parametrization [92], calculated to NLO 
in the MS scheme at a scale Qo = 1.0 GeV 

x(u - u)(x, Ql) 
x(d — d)(x, Ql) 
xs(x, Ql) = xs(x, Ql) 
x(d - u)(x, Ql) 
x(u + d)(x,Ql) 
xg(x,Ql) 

and xqi(x, Ql) = xql(x, Ql) = for qi — c,b,t and we have related the unpolarized input 
distribution to the longitudinally polarized ones by the relations [73] 



= 1.344a; a501 (l - x) 3 - 689 (l + 6.402x - 873 ) 

= 0.64a; - 501 (l - x) 4 - 247 (l + 2.69a; - 333 ) 

= 0.064rr-°- 143 (l - x) 8 ' 041 (l + 6.112a;) 

= 0.071a; a501 (l-a;) 8 - 041 (l + 30.0a;) 

= 0.255a;- - 143 (l - a;) 8 - 041 (l + 6.112a;) 

= 1.123a;- - 206 (l-x) 4 - 673 (l + 4.269a; L508 ) (5.63) 
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xAu(x, Ql) 
xAu(x, Ql) 
xAd(x, Ql) 
xAd{x,Ql) 
xAs(x, Ql) = xAs(x, Ql) 



xrj u (x)xu(x,Ql) 
cos 9 D (x,Ql 
xr] d (x)xd(x, Ql 



x{u — u) — -x(d — d) 



{x, Ql) + xAu(x, Ql) 



= cos6 D (x,Ql) --x(d-d)(x,Ql) +xAd(x,Ql) 
= xr) s (x)xs(x,Ql) 



(5.64) 



and xAqi(x, Ql) = xAqi(x, Ql) = for qi = c, b, t. 

A so-called "spin dilution factor" as defined in [73], which appears in the equations 
above is given by 



cos 9 D (x, Ql) 



1 + 



2a s (Q 2 )(l-x) 



(5.65) 



3 y/x 

In this second (GGR) model, in regard to the initial conditions for the gluons, we have 
made use of two different options, characterized by a parameter r/ dependent on the 
corresponding option. The first option, that we will denote by GGR1, assumes that 
gluons are moderately polarized 



xAg(x,Ql) 
Vu(x) = rj d (x) 
rj s (x) 



x ■ xg(x,Ql) 
-2.49 + 2.8^ 
-1.67 + 2. lv^, 



(5.66) 



while the second option (GGR2) assumes that gluons are not polarized 

xAg(x,Ql) = 
Vu(x) = rj d (x) = -3.03 + 3.0y/x 
rj s (x) = -2.71 + 2.9^- 



(5.67) 



We have plotted both ratios A T /f + and differences (xf + —xA T f) for various flavours as 
a function of x. For the up quark, while the two models GGR1 and GGR2 are practically 
overlapping, the difference between the GGR and the GRSV models in the the ratio 
Atu/u + is only slightly remarked in the intermediate x region (0.1 — 0.5). In any case, 
it is just at the few percent level (Fig. (5.6)), while the inequality is satisfied with 
a ratio between the plus helicity distribution and transverse around 10 percent from 
the saturation value, and above. There is a wider gap in the inequality at small x, 
region characterized by larger transverse distribution, with values up to 40 percent from 
saturation. A similar trend is noticed for the x-behaviour of the inequality in the case 
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of the down quark (Fig. 5.7). In this latter case the GGR and the GRSV model show 
a more remarked difference, especially for intermediate x-values. An interesting features 
appears in the corresponding plot for the strange quark (Fig. (5. 8)), showing a much wider 
gap in the inequality (50 percent and higher) compared to the other quarks. Here we have 
plotted results for the two GGR models (GGR1 and GGR2). Differently from the case 
of the other quarks, in this case we observe a wider gap between lhs and rhs at larger x 
values, increasing as x — > 1. In figs. (5.9)and (5.10) we plot the differences (xf + —xA T f) 
for strange and charm and for bottom and top quarks respectively, which show a much 
more reduced evolution from the saturation value up to the final corresponding evolving 
scales (100 and 200 GeV). As a final application we finally discuss the behaviour of the 
tensor charge of the up quark for the two models as a function of the final evolution scale 
Q. We recall that like the isoscalar and the isovector axial vector charges defined from 
the forward matrix element of the nucleon, the nucleon tensor charge is defined from the 
matrix element of the tensor current 

{PS T \i>G^ l5 \ a i>\P, S T ) = 25q a (Ql) (P»S T - P»S») (5.68) 

where 5 a q(Ql) denotes the flavour (a) contribution to the nucleon tensor charge at a scale 
Qo and St is the transverse spin. 

In fig. (5.11) we plot the evolution of the tensor charge for the models we have taken in 
exam. At the lowest evolution scales the charge is, in these models, above 1 and decreases 
slightly as the factorization scale Q increases. We have performed an evolution up to 200 
GeV as an illustration of this behaviour. There are substantial differences between these 
models, as one can easily observe, which are around 20 percent. From the analysis of 
these differences at various factorization scales we can connect low energy dynamics to 
observables at higher energy, thereby distinguishing between the various models. Inclusion 
of the correct evolution, up to subleading order is, in general, essential. 

5.8 Conclusions 

We have illustrated the use of x-space based algorithms for the solution of evolution 
equations in the leading and in the next-to-leading approximation and we have provided 
some applications of the method both in the analysis of Soffer's inequality and in the 
investigation of other relations, such as the evolution of the proton tensor charge, for 
various models. The evolution has been implemented using a suitable base, relevant for 
an analysis of positivity in LO, using kinetic arguments. The same kinetic argument 
has been used to prove the positivity of the evolution of hi and of the tensor charge 
up to NLO. In our implementations we have completely relied on recursion relations 
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without any reference to Mellin moments. We have provided several illustrations of the 
recursive algorithm and extended it to the nonforward evolution up to NLO. Building on 
previous work for the forward evolution, we have presented a master-form of the nonsinglet 
evolution of the skewed distributions, a simple proof of positivity and a related Kramers 
Moyal expansion, valid in the DGLAP region of the skewed evolution for any value of 
the asymmetry parameter (. We hope to return with a complete study of the nonforward 
evolution and related issues not discussed here in the near future. 



pa 
+ 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 5.4: Coefficients A n (x) + a s (Q 2 )B n , with n — 0, . . . , 4 for a final scale Q 
GeV for the quark up. 
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Figure 5.5: Coefficients A n (x) + a s (Q 2 )B n , with n — 0, . . . , 4 for a final scale Q = 100 
GeV for the quark down. 
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Figure 5.6: Test of Soffer's inequality for quark up at Q = 100 GeV for different models. 
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Figure 5.9: Soffer's inequality for strange and charm in the GRSV model. 
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Figure 5.10: Soffer's inequality for bottom and top in the GRSV model. 
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Figure 5.11: Tensor charge as a function of Q for up and down quark for the GRSV 
and GGR models. 



Chapter 6 



Simulation of air showers: lateral 
distributions, multiplicities and the 
search for new physics at Auger 

6.1 A brief overview 

The inclusive high energy cosmic rays spectrum has been the focus point of a lot of 
attention in the last few years. The spectrum, at a first look, appears to be deprived of 
any structure and, on the other hand, drops dramatically as we move upward in energy. 
The region which we are interested in, in our study, is the high energy region, where a 
correct description of the underlying quark-gluon dynamics is important in order to gain 
some insight into the physics of the extensive air showers that form when a high energy 
primary (say a proton or a neutrino) collides with the nuclei in the atmosphere. Here we 
would like to put things into place in order to set up the stage for our analysis, which will 
be developed in this chapter and will be specialized in the next chapter to the case of mini 
black hole production in cosmic rays in the context of brane models. Our journey has 
been guided by our interest in applying standard QCD tools, such as the RG evolution 
of fragmentation functions, to this challenging environment. It turns out, as expected, 
that the issues to be analyzed, using the parton model, are manifold and do not have 
a clear-cut answer at the moment, although there is a general agreement that models 
based on reggeon theory, popular in the sixties, can capture the complex dynamics of a 
proton-air collision at those energies. The upper end of the cosmic rays spectrum (these 
cosmic rays are classified as ultra high energy cosmic rays or UHECR), in particular, has 
been affected by a controversy which has spanned several decades over the existence of a 
cutoff which would more or less sharply limit it to stay below 10 19 eV. 

In this chapter we perform large scale air shower simulations around the GZK cutoff 
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are performed. An extensive analysis of the behaviour of the various subcomponents 
of the cascade is presented. We focus our investigation both on the study of total and 
partial multiplicities along the entire atmosphere and on the geometrical structure of the 
various cascades, in particular on the lateral distributions. The possibility of detecting new 
physics in Ultra High Energy Cosmic Rays (UHECR) at Auger is also investigated. We 
try to disentangle effects due to standard statistical fluctuations in the first proton impact 
in the shower formation from the underlying interaction and comment on these points. 
We argue that theoretical models predicting large missing energy may have a chance to 
be identified, once the calibration errors in the energy measurements are resolved by the 
experimental collaborations, in measurements of inclusive multiplicities. 

6.2 Introduction 

One of the most intriguing experimental observations of recent years is the detection 
of Ultra-High-Energy-Cosmic-Rays (UHECR), with energy in excess of the Greisen- 
Zatsepin-Kuzmin (GZK) cutoff [77] (for a review see [8]). While its validity is still under 
some dispute, it is anticipated that the forthcoming Auger [11] and EUSO [57] experiments 
will provide enough statistics to resolve the debate. From a theoretical perspective, the 
Standard Model of particle physics and its Grand Unified extensions indicate that many 
physical structures may lie far beyond the reach of terrestrial collider experiments. If this 
eventuality materializes it may well be that the only means of unlocking the secrets of the 
observed world will be mathematical rigor and peeks into the cosmos in its most extreme 
conditions. In this context the observation of UHECR is especially puzzling because of 
the difficulty in explaining the events without invoking some new physics. There are 
apparently no astrophysical sources in the local neighborhood that can account for the 
events. The shower profile of the highest energy events is consistent with identification of 
the primary particle as a hadron but not as a photon or a neutrino. The ultrahigh energy 
events observed in the air shower arrays have muonic composition indicative of hadrons. 
The problem, however, is that the propagation of hadrons over astrophysical distances is 
affected by the existence of the cosmic background radiation, resulting in the GZK cutoff 
on the maximum energy of cosmic ray nucleons -Egzk < 10 20 eV [76]. Similarly, photons 
of such high energies have a mean free path of less than 10 Mpc due to scattering from the 
cosmic background radiation and radio photons. Thus, unless the primary is a neutrino, 
the sources must be nearby. On the other hand, the primary cannot be a neutrino because 
the neutrino interacts very weakly in the atmosphere. A neutrino primary would imply 
that the depths of first scattering would be uniformly distributed in column density, which 
is contrary to the observations. 
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The most exciting aspect of the UHECR is the fact that the Auger and EUSO ex- 
periments will explore the physics associated with these events, and provide a wealth of 
observational data. Clearly, the first task of these experiments is to establish whether 
the GZK cutoff is violated, and to settle the controversy in regard to the air shower 
measurement. 

6.3 Probing new physics with UHECR 

We may, however, entertain the possibility that these experiments can probe various 
physics scenarios. In the first place, the center of mass energy in the collision of the 
primary with the atmosphere is of the order of 100 TeV and exceed the contemporary, 
and forthcoming, collider reach by two orders of magnitude. Thus, in principle the air 
shower analysis should be sensitive to any new physics that is assumed to exist between 
the electroweak scale and the collision scale due to the interaction of the primaries with 
the atmospheric nuclei. Other exciting possibilities include the various explanations that 
have been put forward to explain the existence of UHECR events [8, 112, 17, 16, 18, 43], 
and typically assume some form of new physics. One of the most intriguing possible 
solutions is that the UHECR primaries originate from the decay of long-lived super- 
heavy relics, with mass of the order of 10 12 ~ 15 GeV [17, 16, 18, 43]. In this case the 
primaries for the observed UHECR would originate from decays in our galactic halo, and 
the GZK bound would not apply. This scenario is particularly interesting due to the 
possible connection with superstring theory. From the particle physics perspective the 
meta-stable super-heavy candidates should possess several properties. First, there should 
exist a stabilization mechanism which produces the super-heavy state with a lifetime of 
the order of 10 17 s < r x < 10 28 s, and still allows it to decay and account for the observed 
UHECR events. Second, the required mass scale of the meta-stable state should be of 
order, Mx ~ 10 12 ~ 13 GeV. Finally, the abundance of the super-heavy relic should satisfy 
the relation (flx/^o){to/ T x) ~ 5 x 10~ n , to account for the observed flux of UHECR 
events. Here t is the age of the universe, Tx the lifetime of the meta-stable state, Q is 
the critical mass density and fix is the relic mass density of the meta-stable state. It is 
evident that the parameters of the super-heavy meta-stable states are sufficiently flexible 
to accommodate the observed flux of UHECR, while evading other constraints [43]. 

Superstring theory inherently possesses the ingredients that naturally give rise to 
super-heavy meta-stable states. Such states arise in string theory due to the breaking 
of the non-Abelian gauge symmetries by Wilson lines. The massless spectrum then con- 
tains states with fractional electric charge or "fractional" U(l)z> charge [121, 54, 35]. The 
lightest states are meta-stable due to a local gauge, or discrete, symmetry [35, 58]. This 
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phenomenon is of primary importance for superstring phenomenology. The main conse- 
quence is that it generically results in super-massive states that are meta-stable. The 
super-heavy states can then decay via the nonrenormalizable operators, which are pro- 
duced from exchange of heavy string modes, with lifetime r x > 10 7 ~ 17 years [54, 16, 43]. 
The typical mass scale of the exotic states will exceed the energy range accessible to 
future collider experiments by several orders of magnitude. The exotic states are ren- 
dered super-massive by unsuppressed mass terms [59], or are confined by a hidden sector 
gauge group [54]. String models may naturally produce mass scales of the required or- 
der, Mx ~ 10 12_13 GeV, due to the existence of an hidden sector that typically contains 
non-Abelian SU{n) or SO(2n) group factors. The hidden sector dynamics are set by 
the initial conditions at the Planck scale, and by the hidden sector gauge and matter 
content, M x ~ A^-^(N,n f ). Finally, the fact that M x ~ 10 12 ~ 13 GeV implies that 
the super-heavy relic is not produced in thermal equilibrium and some other production 
mechanism is responsible for generating the abundance of super-heavy relic [35]. 

The forthcoming cosmic rays observatories can therefore provide fascinating experi- 
mental probes, both to the physics above the electroweak scale as well as to more exotic 
possibilities at a much higher scale. It is therefore imperative to develop the theoretical 
tools to decipher the data from these experiments. Moreover, improved information on 
the colliding primaries may reveal important clues on the properties of the decaying meta- 
stable state, which further motivates the development of such techniques. In this chapter 
we make a modest step in this direction, by studying possible modifications of air shower 
simulations, that incorporate the possible effects of new physics above the electroweak 
scale. This is done by varying the cross section in the air shower codes that are used by 
the experimentalists. In this respect we assume here for concreteness that the new physics 
above the electroweak scale remains perturbative and preserve unitarity, as in the case 
of supersymmetric extensions of the Standard Model. This in turn is motivated by the 
success of supersymmetric gauge coupling unification [55] and their natural incorporation 
in string theories. In the case of supersymmetry the deviations from the Standard Model 
are typically in the range of a few percent, a quantitative indication which we take as our 
reference point for study. 

6.3.1 Possible Developments 

Even if the forthcoming experiments will confirm the existence of UHECR events, it 
remains to be seen whether any new physics can be inferred from the results. We will 
argue that this is a very difficult question. 

A possible way, in the top-down models of the UHECR interaction is to optimize 
the analysis of any new high energy primary interaction. One should keep in mind that 
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the information carried by the primaries in these collisions is strongly "diluted" by their 
interaction with the atmosphere and that large statistical fluctuations are immediately 
generated both by the randomness of the first impact, the variability in the zenith angle of 
the impact, and the natural fluctuations in the - extremely large - phase space available at 
those energies. We are indeed dealing with extreme events. These uncertainties are clearly 
mirrored even in the existing Monte Carlo codes for the simulation of air showers, and, 
of course, in the real physical process that these complex Monte Carlo implementations 
try, at their best, to model (see also [6] for the discussion of simulation issues) . Part of 
our work will be concerned primarily with trying to assess, by going through extensive 
air-shower simulations using existing interaction models - at the GZK and comparable 
energies - the main features of the showers, such as the multiplicities at various heights 
and on the detector plane. We will illustrate the geometry of a typical experimental setup 
to clarify our method of analysis and investigate in detail some geometrical observables. 

A second part of our analysis will be centered around the implications of a modified 
first impact on the multiplicities of the subcomponents. Our analysis here is just a first 
step in trying to see whether a modified first impact cross section has any implication on 
the multiplicity structure of the shower. The analysis is computationally very expensive 
and has been carried out using a rather simple strategy to render it possible. We critically 
comment on our results, and suggest some possible improvements for future studies. 

6.4 Simulation of air showers 

The quantification of the variability and parametric dependence of the first impact in the 
formation of extensive air showers can be discussed, at the moment, only using Monte 
Carlo event generators. Although various attempts have been made in the previous liter- 
ature to model the spectrum of a generic "X-particle" decay in various approximations, all 
of them include - at some level and with variants - some new physics in the generation of 
the original spectrum. In practice what is seen at experimental level is just a single event, 
initiated by a single hadron (a proton) colliding with an air nucleus (mostly of oxygen or 
nitrogen) within the 130 km depth of the Earth atmosphere. Our studies will show that 
the typical strength of the interaction of the primary at the beginning of the showers - at 
least using the existing Monte Carlo codes - has to increase fairly dramatically in order 
to be able to see - at the experimental level - any new physics. 

Our objective here is to assess the actual possibility, if any, to detect new physics from 
the high energy impact of the primary cosmic ray assuming that other channels open up 
at those energies. Our investigation here is focused on the case of supersymmetry, which 
is the more widely accepted extension of the Standard Model. Other scenarios are left for 
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future studies. 

We recall that at the order of the GZK cutoff, the center of mass energy of the first 
collision reaches several hundreds of TeVs and is, therefore, above any supersymmetric 
scale, according to current MSSM models. It is therefore reasonable to ask weather 
supersymmetric interactions are going to have any impact on some of the observables 
that are going to be measured. 

We will provide enough evidence that supersymmetric effects in total hadronic cross 
sections cannot raise the hadronic nucleon nucleon cross section above a (nominal) 100% 
upper limit. We will then show that up to such limit the fluctuations in 1) the multiplicity 
distributions of the most important components of the (ground) detected air showers and 
2) the geometric distributions of particles on the detector are overwhelmingly affected by 
natural (statistical) fluctuations in the formation of the air showers and insignificantly by 
any interaction whose strength lays below such 100% nominal limit. 

In order to proceed with our analysis we need to define a set of basic observables which 
can be used in the characterization of the shower at various heights in the atmosphere. 

There are some basic features of the shower that are important in order to understand 
its structure and can be summarized in: 1) measurements of its multiplicities in the main 
components; 2) measurements of the geometry of the shower. Of course there are obvious 
limitations in the study of the development of the shower at the various levels, since the 
main observations are carried out on the ground. However, using both Cerenkov telescopes 
and fluorescence measurements by satellites one hopes to reconstruct the actual shape of 
the shower as it develops in the atmosphere. 

To illustrate the procedure that we have implemented in order to characterize the 
shower, we have assumed that the first (random) impact of the incoming primary (proton) 
cosmic ray takes place at zero zenith angle, for simplicity. We have not carried out 
simulations at variable zenith, since our objective is to describe the main features of the 
shower in a rather simple, but realistic, setting. We have chosen a flat model for the 
atmosphere and variable first impacts, at energies mainly around the GZK cutoff region. 
Our analysis has been based on CORSIKA [78] and the hadronization model chosen has 
been QGSJET [83]. 

Measurements at any level are performed taking the arrival axis (z-axis) of the primary 
as center of the detector. The geometry of the shower on the ground and at the various 
selected observation levels has been always measured with respect to this axis. The 
"center" of the detector is, in our simulations, assumed to be the point at which the z-axis 
intersects the detector plane. Below, the word "center" refers to this particular geometrical 
setting. 

The shower develops according to an obvious cylindrical symmetry around the vertical 
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z-axis, near the center. The various components of the showers are characterized at any 
observation level by this cylindrical symmetry. Multiplicities are plotted after integration 
over the azimuthal angle and shown as a function of the distance from the core (center), 
in the sense specified above. 

The showers show for each subcomponent specific locations of the maxima and widths 
of the associated distributions. We will plot the positions of the maxima along the entire 
spatial extent of the shower in the atmosphere. These plots are useful in order to have an 
idea of what is the geometry of the shower in the 130 km along which it develops. 

6.5 Features of the Simulation 

Most of our simulations are carried out at two main energies, 10 19 and 10 20 eV. Simulations 
have been performed on a small cluster running a communication protocol (openmosix) 
which distributes automatically the computational load. The simulation program fol- 
lows each secondary from beginning to end and is extremely time and memory intensive. 
Therefore, in order to render our computation manageable we have implemented in COR- 
SIKA the thinning option [79], which allows to select only a fraction of the entire shower 
and followed its development from start to end. We recall that CORSIKA is, currently, 
the main program used by the experimental collaborations for the analysis of cosmic rays. 
The results have been corrected statistically in order to reproduce the result of the actual 
(complete) shower. The CORSIKA output has been tokenized and then analyzed using 
various intermediate software written by us. The number of events generated, even with 
the thinning algorithm, is huge at the GZK energy and requires an appropriate handling of 
the final data. We have performed sets of run and binned the data using bins of 80 events, 
where an event is a single impact with its given parametric dependence. The memory 
cost of a statistically significant set of simulations is approximately 700 Gigabytes, having 
selected in our simulations a maximal number of observations levels (9) along the entire 
height of the atmosphere. 

6.5.1 Multiplicities on the Ground 

We show in Fig. 2.36 results for the multiplicities of the photon component and of the e ± 
components at 10 19 eV plotted against the distance from the core (center) of the detector. 
For photons, the maximum of the shower is around 90 meters from the center, as measured 
on the plane of the detector. As evident from the plot, the statistics is lower as we get 
closer to the central axis (within the first 10 meters from the vertical axis), a feature 
which is typical of all these distributions, given the low multiplicities measured at small 
distances from the center. For electrons and positrons the maxima also lie within the first 
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100 meters, but slightly closer to the center and are down by a factor of 10 in multiplicities 
with respect to the photons. Positron distributions are suppressed compared to electron 
distributions. It is also easily noticed that the lower tail of the photon distribution is 
larger compared to the muon distribution, but all the distributions show overall similar 
widths, about 1 km wide. 

Multiplicities for muons (see Fig. 2.37) are a factor 1000 down with respect to photons 
and 100 down with respect to electrons. The maxima of the muon distributions are also 
at comparable distance as for the photons, and both muon and antimuons show the same 
multiplicity. 

At the GZK energy (Figs 2.38, 2.39) the characteristics of the distributions of the 
three main components (photons, electrons, muons) do not seem to vary appreciably, 
except for the values of the multiplicities, all increased by a factor of 10 respect to the 
previous plots. The maxima of the photons multiplicities are pushed away from the center, 
together with their tails. There appears also to be an increased separation in the size of 
the multiplicities of electrons and positrons and a slightly smaller width for the photon 
distribution compare to the lower energy result (10 19 eV). We should mention that all 
these gross features of the showers can possibly be tested after a long run time of the 
experiment. Our distributions have been obtained averaging over sets of 80 events with 
independent first impacts. 

6.6 Missing multiplicities? 

Other inclusive observables which are worth studying are the total multiplicities, as mea- 
sured at ground level, versus the total energy of the primary. We show in Fig. 6.5 a double 
logarithmic plot of the total multiplicities of the various components versus the primary 
energy in the range 10 15 — 10 20 eV, which appears to be strikingly linear. From our result 
it appears that the multiplicities can be fitted by a relation of the form y=m*x+q, where 
y = Log(N) and x = Log(E) or N(E) = 10 q x m . The values of m and q are given by 
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where m is the slope, m appears to be almost universal for all the components, while 
the intercept (q) depends on the component. The photon component is clearly dominant, 
followed by the electron, positron and the two muon components which appear to be 
superimposed. It would be interesting to see whether missing energy effects, due, for 
instance, to an increased multiplicity rate toward the production of weakly interacting 
particles can modify this type of inclusive measurements, thereby predicting variations in 
the slopes of the multiplicities respect to the Monte Carlo predictions reported here. One 
could entertain the possibility that a failure to reproduce this linear behaviour could be a 
serious problem for the theory and a possible signal of new physics. Given the large sets of 
simulations that we have performed, the statistical errors on the Monte Carlo results are 
quite small, and the Monte Carlo prediction appear to be rather robust. The difficulty of 
these measurements however, lay mainly in the energy reconstruction of the primary, with 
the possibility of a systematic error. However, once the reconstruction of the energy of the 
primary is under control among the various UHECR detectors, with a global calibration, 
these measurements could be a possible test for new physics. At the moment, however, 
we still do not have a quantification of the deviation from this behaviour such as that 
induced by supersymmetry or other competing theoretical models. 

6.6.1 Directionality of the Bulk of the Shower 

Another geometrical feature of the shower is the position of its bulk (measured as the 
opening of the radial cone at the radial distance where the maximum is achieved) as 
a function of the energy. This feature is illustrated in Fig. 6.6. Here we have plotted 
the averaged location of the multiplicity distribution as a function of the energy of the 
incoming primary. The geometrical center of the distributions tend to move slightly 
toward the vertical axis (higher collimation) as the energy increases. From the same plot 
it appears that the distributions of electrons, positrons and photons are closer to the 
center of the detector compared to the muon-antimuon distributions. As shown in the 
figure, the statistical errors on these results appear to be rather small. 

6.6.2 The Overall Geometry of the Shower 

As the shower develops in the atmosphere, we can monitor both the multiplicities of the 
various components and the average location of the bulk of the distributions at various 
observation levels. As we have mentioned above, we choose up to 9 observation levels 
spaces at about 13 km from one another. The lowest observation level is, according to 
our conventions, taken to coincide with the plane of the detector. 

We show in Fig. 6.7a a complete simulation of the shower using a set of 9 observation 
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levels, as explained above, at an energy of 10 eV. In the simulation we assume that a 
first impact occurs near the top of the atmosphere at an eight of 113 km and we have 
kept this first impact fixed. The multiplicities show for all the components a rather fast 
growth within the first 10 km of crossing of the atmosphere after the impact, with the 
photon components growing faster compared to the others. The electron component also 
grows rather fast, and a similar behaviour is noticed for the muon/antimuon components, 
which show a linear growth in a logarithmic scale (power growth). In the following 40 
km downward, from a height of 100 km down to 60 km, all the components largely 
conserve their multiplicities. Processes of regeneration of the various components and 
their absorption seem to balance. For the next 20 km, from a height of 60 km down to 40 
km, all the components starts to grow, with the photon component showing a faster (power 
growth) with the traversed altitude. Slightly below 40 km of altitude the multiplicities of 
three components seem to merge (muons and electrons), while the photon component is 
still dominant by a factor of 10 compared to the others. The final development of the air 
shower is characterized by a drastic growth of all the components, with a final reduced 
muon component, a larger electron component and a dominant photon component. The 
growth in this last region (20 km wide) and in the first 20 km after the impact of the 
primary -in the upper part of the atmosphere- appear to be comparable. The fluctuations 
in the multiplicities of the components are rather small at all levels, as shown (for the 
photon case) in Fig. 6. 7b. 



As we reach the GZK cutoff, increasing the energy of the primary by a factor of 10, 
the pattern just discussed in Fig. 6.7a is reproposed in Fig. 6.8a, though -in this case- the 
growth of the multiplicities of the subcomponents in the first 20 km from the impact and 
in the last 20 km is much stronger. The electron and the muon components appear to be 
widely separated, while the electron and positron components tend to be more overlapped. 
To the region of the first impact -and subsequent growth- follows an intermediate region, 
exactly as in the previous plots, where the two phenomena of production and absorption 
approximately balance one another and the multiplicities undergo minor variations. The 
final growth of all the components is, a this energy, slightly anticipated compared to Fig. 
6.7a, and starts to take place at a height of 40 km and above and continues steadily until 
the first observation level. The photon remains the dominant component, followed by the 
electron and the muon component. Also in this case the fluctuations (pictured for the 
photon component only, Fig. 6.8b) are rather small. 
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6.6.3 The Opening of the Shower 

In our numerical study the geometrical center of each component of the shower is identified 
through a simple average with respect to all the distances from the core 



where N is the total multiplicity at each selected observation level, Ri is the position of 
the produced particle along the shower and i runs over the single events. This analysis 
has been carried out for 9 equally spaced levels and the result of this study are shown in 
Figs. 6.9,6.10 and 6.8 at two different values of energy (10 19 and 10 20 eV). The opening 
of the various components are clearly identified by these plots. We start from Fig. 6.9. 
We have taken in this figure an original point of impact at a height of 113 Km, as in 
the previous simulations. It is evident that the photon component of the shower tends 
to spread rather far and within the first 20 km of depth into the atmosphere has already 
reached an extension of about 2 km; reaching a lateral extension of 10 km within the first 
60 km of crossing of the atmosphere. 

Starting from a height of 50 km down to 10 km, the shower gets reabsorbed (turns 
toward the center) and is characterized by a final impact which lays rather close to the 
vertical axis. Electrons and photons follow a similar behaviour, except that for electrons 
whose lateral distribution in the first section of the development of the shower is more 
reduced. The muonic (antimuons) subcomponents appear to have a rather small opening 
and develop mostly along the vertical axis of impact. In the last section of the shower all 
the components get aligned near the vertical axis and hit the detector within 1 km. 

Few words should be said about the fluctuations. At 10 19 eV, as shown in Fig. 6.9b, the 
fluctuations are rather large, especially in the first part of the development of the shower. 
These turn out to be more pronounced for photons, whose multiplicity growth is large and 
very broad. As we increase the energy of the primary to 10 20 eV the fluctuations in the 
lateral distributions (see Fig. 6.10b) are overall reduced, while the lateral distributions of 
the photons appear to be drastically reduced (Fig. 6.10a). 

6.7 Can we detect new physics at Auger? 

There are various issues that can be addressed, both at theoretical and at experimental 
level, on this point, one of them being an eventual confirmation of the real existence of 
events above the cutoff. However, even if these measurement will confirm their existence, 
it remains yet to be seen whether any additional new physics can be inferred just from an 
analysis of the air shower. A possibility might be supersymmetry or any new underlying 
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interaction, given the large energy available in the first impact. We recall that the spec- 
trum of the decaying X-particle (whatever its origin may be), prior to the atmospheric 
impact of the UHECR is of secondary relevance, since the impact is always due to a single 
proton. Unless correlations are found among different events - and by this we mean that 
a large number of events should be initiated by special types of primaries - we tend to 
believe that effects due to new interactions are likely to play a minor role. 

In previous works we have analyzed in great detail the effects of supersymmetry in 
the formation of the hadronic showers. These studies, from our previous experience [40, 
41, 42], appear to be rather complex since they involve several possible intermediate and 
large final scales and cannot possibly be conclusive. There are some obvious doubts that 
can be raised over these analysis, especially when the DGLAP equations come to be 
extrapolated to such large evolution scales, even with a partial resummation of the small- 
x logarithms. In many cases results obtained in this area of research by extrapolating 
results from collider phenomenology to extremely high energies should be taken with 
extreme caution in order to reduce the chances of inappropriate hasty conclusions. What 
is generally true in a first approximation is that supersymmetric effects do appear to 
be mild [40, 41, 42]. Rearrangements in the fragmentation spectra or supersymmetric 
effects in initial state scaling violations are down at the few percent level. We should 
mention that the generation of supersymmetric scaling violations in parton distributions, 
here considered to be the bulk of the supersymmetric contributions, are rather mild if the 
entrance into the SUSY region takes place "radiatively" as first proposed in [40, 41]. This 
last picture might change in favour of a more substantial signal if threshold enhancements 
are also included in the evolution, however this and other related points have not yet been 
analyzed in the current literature. 

6.7.1 The Primary Impact and a Simple Test 

Our objective, at this point, is to describe the structural properties of the shower with an 
emphasis on the dynamics of the first impact of the primary with the atmosphere, and 
at this stage one may decide to look for the emergence of possible new interactions, the 
most popular one being supersymmetry. 

One important point to keep into consideration is that the new physical signal carried 
by the primaries in these collisions is strongly "diluted" by their interaction with the 
atmosphere and that large statistical fluctuations are immediately generated both by the 
randomness of the first impact, the variability in the zenith angle of the impact, and 
the -extremely large- phase space available at those energies in terms of fragmentation 
channels. We can't possibly underestimate these aspects of the dynamics, which are 
at variance with previous analysis, where the search for supersymmetric effects (in the 
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vacuum) seemed to ignore the fact that our detectors are on the ground and not in space. 

For this reason we have resorted to a simple and realistic analysis of the structure of 
air showers as can be obtained from the current Monte Carlo. 

The simplest way to test whether a new interaction at the first proton-proton impact 
can have any effect on the shower is to modify the cross section at the first atmospheric 
impact using CORSIKA in combination with some of the current hadronization models 
which are supposed to work at and around the GZK cutoff. There are obvious limitations 
in this approach, since none of the existing codes incorporates any new physics beyond 
the standard model, but this is possibly one of the simplest ways to proceed. For this 
purpose we have used SYBILL [62], with the appropriate modifications discussed below. 

To begin let's start by recalling one feature of the behaviour of the hadronic cross 
section (pp or pp) at asymptotically large energies. There is evidence (see [61]) demon- 
strating a saturation of the Froissart bound of the total cross section with rising total 
energy, s. This log 2 (s) growth of the total cross section is usually embodied into many of 
the hadronization models used in the analysis of first impact and leaves, therefore, little 
room for other substantial growth with the opening of new channels, supersymmetry be- 
ing one of them. We should also mention that various significant elaborations [19] on the 
growth of the total cross section and the soft pomeron dominance have been discussed 
in the last few years and the relation of this matter [123] with the UHECR events is of 
utmost relevance. 

With this input in mind we can safely "correct" the total cross section by at most a 
(nominal) factor of 2 and study whether these nominal changes can have any impact on 
the structural properties of the showers. 

We run simulations on the showers generated by this modification and try to see 
whether there is any signal in the multiplicities which points toward a structural (multi- 
plicity, geometrical) modification of the air showers in all or some of its subcomponents. 
For this purpose we have performed runs at two different energies, at the GZK cutoff and 
1 decade below, and analyzed the effects due to these changes. 

We show in two figures results on the multiplicities, obtained at zero zenith angle, of 
some selected particles (electrons and positrons, in our case, but similar results hold for 
all the dominant components of the final shower) obtained from a large scale simulation of 
air showers at and around the GZK cutoff. We have used the simulation code CORSIKA 
for this purpose. 

In Figs. 6.11-6.15 we show plots obtained simulating an artificial first proton impact 
in which we have modified the first interaction cross section by a nominal factor ranging 
from 0.7 to 2. We plot on the y-axis the corresponding fluctuations in the multiplicities 
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both for electrons and positrons. Statistical fluctuations 1 have been estimated using bins 
of 80 runs. The so-developed showers have been thinned using the Hillas algorithm, as 
usually done in order to make the results of these simulations manageable, given the size 
of the showers at those energies. As one can immediately see, the artificial corrections on 
the cross section are compatible with ordinary fluctuations of the air-shower. We have 
analyzed all the major subcomponents of the air shower, photons and leptons, together 
with the corresponding neutrino components. We can summarize these findings by as- 
sessing that a modified first impact, at least for such correction factors in the 0.7-2 range 
in the cross section, are unlikely to modify the multiplicities in any appreciable way. 

A second test is illustrated in Figs. 6.16-6.20. Here we plot the same correction factors 
on the x-axis as in the previous plots (6.11-6.15) but we show on the y-axis (for the same 
particles) the average point of impact on the detector and its corresponding statistical 
fluctuations. As we increase the correction factor statistical fluctuations in the formation 
of the air shower seem to be compatible with the modifications induced by the "new 
physics" of the first impact and no special new effect is observed. 

Fluctuations of these type, generated by a minimal modification of the existing codes 
only at the first impact may look simplistic, and can possibly be equivalent to ordinary 
simulations with a simple rescaling of the atmospheric height at which the first collision 
occurs, since the remaining interactions are, in our approach, unmodified. The effects 
we have been looking for, therefore, appear subleading compared to other standard fluc- 
tuations which take place in the formation of the cascade. On the other hand, drastic 
changes in the structure of the air shower should possibly depend mostly on the physics 
of the first impact and only in a less relevant way on the modifications affecting the cas- 
cade that follows up. We have chosen to work at an energy of 10 20 eV but we do not 
observe any substantial modifications of our results at lower energies (10 19 eV), except for 
the multiplicities which are down by a factor of 10. Our brief analysis, though simple, 
has the purpose to illustrate one of the many issues which we believe should be analyzed 
with great care in the near future: the physics of the first impact and substantial addi- 
tional modifications to the existing codes in order to see whether any new physics can be 
extracted from these measurements. 

6.8 Conclusions 

We have tried to analyze with a searching criticism the possibility of detecting new physics 
at Auger using current ideas about supersymmetry, the QCD evolution and such. While 

1 we keep the height of the first proton impact with the atmosphere arbitrary for each selected correction 
factor (x-axis) 
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the physics possibilities of these experiments are far reaching and may point toward a 
validation or refusal of the existence of a GZK cutoff, we have argued that considerable 
progress still needs to be done in order to understand better the hadronization models at 
very large energy scales. Our rather conservative viewpoint stems from the fact that the 
knowledge of the structure of the hadronic showers at large energies is still under debate 
and cannot be conclusive. We have illustrated by an extended and updated simulation 
some of the characteristics of the showers, the intrinsic fluctuations in the lateral distri- 
butions, the multiplicities of the various subcomponents and of the total spectrum, under 
some realistic conditions. We have also tried to see whether nominal and realistic changes 
in the cross section of the first impact may affect the multiplicities, with a negative out- 
come. We have however pointed out, in a positive way, that new physics models predicting 
large missing energy may have a chance to be identified, since the trend followed by the 
total multiplicities (in a log-log scale) appear to be strikingly linear. 

We do believe, however, that other and even more extensive simulation studies should 
be done, in combination with our improved understanding of small-x effects in QCD 
at large parton densities, in order to further enhance the physics capabilities of these 
experiments. Another improvement in the extraction of new physics signals from the 
UHECR experiments will come from the incorporation of new physics in the hadronization 
models. 
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Figure 6.1: Multiplicities of photons and e at the ground level with a proton primary of 
10 19 eV as a function of the distance from the core of the shower. 
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Figure 6.2: Multiplicities of (i at the ground level with a proton primary of 10 19 eV as a 
function of the distance from the core of the shower. 
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Figure 6.3: Multiplicities of photons and e at the ground level with a proton primary of 
10 20 eV as a function of the distance from the core of the shower. 
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Figure 6.4: Multiplicities of (i at the ground level with a proton primary of 10 20 eV as a 
function of the distance from the core of the shower. 
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Figure 6.6: Average distance from the core of the shower of photons, e , n at the ground 
level function of the primary energy. 
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Figure 6.7: Multiplicities of photons, e , n at various levels of observations for a primary 
energy of 10 19 eV. The first impact is forced to occur at the top of the atmosphere. In 
the subfigure (b) we show the uncertainties just in the case of the photons. 
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Figure 6.8: Multiplicities of photons, e ± , ^ at various levels of observations for a primary 
energy of 10 20 eV. The first impact is forced to occur at the top of the atmosphere. In 
the subfigure (b) we show the uncertainties just in the case of the photons. 
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Figure 6.9: Average core distances of photons, e , ji at various levels of observations for a 
primary energy of 10 19 eV. The first impact is forced to occur at the top of the atmosphere. 
In the subfigure (b) we show the uncertainties just in the case of the photons. 
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Figure 6.10: Average core distances of photons, e , ji at various levels of observations 
for a primary energy of 10 20 eV. The first impact is forced to occur at the top of the 
atmosphere. In the subfigure (b) we show the uncertainties just in the case of the photons. 
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Figure 6.11: variation of the photon multiplicity as a function of the first impact cross 
section at 10 19 and at 10 20 eV 
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Figure 6.12: variation of the e ± multiplicity as a function of the first impact cross section 
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Figure 6.13: variation of the // ± multiplicity as a function of the first impact cross section 
at 10 19 and at 10 20 eV 
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Figure 6.15: variation of the multiplicity as a function of the first impact cross section 
at 10 19 and at 10 20 eV 
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Figure 6.16: Lateral distributions of photons as a function of the first impact cross section 
at 10 19 and at 10 20 eV 
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Figure 6.17: Lateral distributions of e ± as a function of the first impact cross section at 
10 19 and at 10 20 eV 



180 



6.8. Conclusions 



19 



Primary energy : 10 eV 



38800 



I 38400 




o 

60 

s- 
O 



37200 



0.7 0.8 0.9 



1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 
cross section (artificial/theoretical) 

20 

Primary energy: 10 eV 



37000 



36800 



36600 



36400 



36200 



36000 



35800 



35600 




0.7 0.8 0.9 



1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1. 

cross section (artificial/theoretical) 



Figure 6.18: Lateral distributions of ^ as a function of the first impact cross section at 
10 19 and at 10 20 eV 
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Figure 6.19: Lateral distributions of v e as a function of the first impact cross section at 
10 19 and at 10 20 eV 
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Figure 6.20: Lateral distributions of as a function of the first impact cross section at 
10 19 and at 10 20 eV 
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Searching for extra dimensions in high 

energy cosmic rays: 

the fragmentation of mini black holes 

7.1 Background on the topic 

Extra Dimensions have been around the theory landscape for some time, starting from the 
20's with the work of Kaluza and Klein. Their original idea was to use a five dimensional 
manifold, the fifth dimension being a circle, and identify the component of the metric 
in this 5-dimensional space as a 4-dimensional vector field. The field, after a Fourier 
expansion in the extra variable, has a massless component which can be identified with 
the photon, plus an infinite tower of massive states. The masses of these additional states 
are inversely proportional to the size of the extra dimensions and their disappearance 
from the low energy spectrum can be simply obtained by tuning the extra dimensions 
appropriately. Then, the natural question to ask is: how large can the extra dimensions 
be without getting into conflict with the experiment? The answer to this question requires 
the Standard Model fields to be localized on the brane with gravity free to propagate in 
the extra space, called "the bulk". This scenario is characterized by a rich phenomenology, 
as we are going to discuss next, with implications which may appear also in the structure 
of cosmic ray showers. One possibility is the formation, due to a lower gravity scale, of 
mini black holes which decay into an "s-wave" or isotropically into all the particles of 
the Standard Model. For this reason we present a study of the main observables of the 
air showers initiated by an incoming primary which collides with the atmosphere and 
characterized by the formation of an intermediate mini black hole. We study particle 
multiplicities, lateral distributions and the ratio of the electromagnetic to the hadronic 
components of these special air showers. In this chapter we illustrate a simulation study 
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of the resulting cascades over the entire range (10 15 — 10 19 eV) of ultra high initial energies, 
for several values of the number of large extra dimensions, for a variety of altitudes of the 
initial interaction and with the energy losses in the bulk taken into account. The results 
are compared with a representative of the standard events, namely the shower due to the 
collision of a primary proton with a nucleon in the atmosphere. Both the multiplicities 
and the lateral distribution of the showers show important differences between the two 
cases and, consequently, may be useful for the observational characterization of the events. 
The electromagnetic/hadronic ratio is strongly fluctuating and, thus, less decisive for the 
altitudes considered. 

This chapter is based on the paper [32]. 

7.2 Introduction 

In the almost structureless fast falling with energy inclusive cosmic ray spectrum, two 
kinematic regions have drawn considerable attention for a long time [8, 112]. These regions 
are the only ones in which the spectral index of the cosmic ray flux shows a sharper 
variation as a function of energy, probably signaling some "new physics", according to 
many. These two regions, termed the knee and the ankle [77] have been puzzling theorists 
and experimentalists alike and no clear and widely accepted explanation of this unusual 
behaviour in the propagation of the primaries - prior to their impact with the earth 
atmosphere - exists yet. A large experimental effort [11, 57] in the next several years will 
hopefully clarify several of the issues related to this behaviour. 

While the ankle is mentioned in the debate regarding the possible existence of the 
so called Greisen, Zatsepin and Kuzmin (GZK) cutoff [76], due to the interaction of the 
primaries with the cosmic background radiation, the proposed resolutions of this puzzle 
are several, ranging from a resonant Z-burst mechanism [120] to string relics and other 
exotic particle decays [17, 16, 18, 43]. The existence of data beyond the cutoff has also 
been critically discussed [13]. 

Given the large energy involved in the first stage of the formation of the air showers, 
the study of the properties of the cascade should be sensitive to any new physics between 
the electroweak scale and the original collision scale. Especially in the highest energy 
region of the spectrum, the energy available in the interaction of the primaries with the 
atmospheric nuclei is far above any conceivable energy scale attainable at future ground- 
based accelerators. Therefore, the possibility of detecting supersymmetry, for instance, in 
cosmic ray showers has also been contemplated [42]. Thus, it is not surprising, that most 
of the attempts to explain these features of the cosmic ray spectrum typically assume 
some form of new physics at those energies. 
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With the advent of theories with a low fundamental scale of gravity [10] and large 
compact or non-compact extra dimensions, the possibility of copiously producing mini 
black holes (based on Thome's hoop conjecture [111]) in collisions involving hadronic 
factorization scales above 1 TeV has received considerable attention [48] [68] and these 
ideas, naturally, have found their way also in the literature of high energy cosmic rays 
[7, 60] and astrophysics [37]. For instance, it was recently suggested that the long known 
Centauro events might be understood as evaporating mini black holes, produced by the 
collision of a very energetic primary (maybe a neutrino) with a nucleon (quark) in the 
atmosphere [99]. Other proposals [69] also either involve new forms of matter (for example 
strangelets) or speculate about major changes in the strong interaction dynamics [124]. 

While estimates for the frequencies of these types of processes both in cosmic rays 
[1, 7, 99] and at colliders [68] have been presented, detailed studies of the multiplicities 
of the particles collected at the detectors, generated by the extensive atmospheric air 
showers following the first impact of the primary rays, are far from covering all the main 
features of the cascade [125]. These studies will be useful in order to eventually disentangle 
new physics starting from an analysis of the geometry of the shower, of the multiplicity 
distributions of its main sub-components [27] and of its directionality from deep space. 
For instance, the study of the location of the maxima of the showers at positions which 
can be detected by fluorescence mirrors [5], generated as they go across the atmosphere, 
and their variations as a function of the parameters of the underlying physical theory, may 
help in this effort [1]; other observables which also contain potential new information are 
the multiplicities of the various particle sub-components and the opening of the showers 
as they are detected on the ground [27]. We will focus on this last type of observables. 

To summarize: in the context of the TeV scale gravity with large extra dimensions 
it is reasonable to assume that mini black holes, black holes with mass of a few TeV, 
can form at the first impact of ultra high energy primary cosmic rays with nucleons in 
the atmosphere. The black hole will evaporate into all types of particles of the Standard 
Model and gravity. The initial partons will hadronize and all resulting particles as they 
propagate in the atmosphere will develop into a shower (s), which eventually will reach the 
detectors. The nature and basic characteristics of these showers is the question that is the 
main subject of the present work. What is the signature on the detector of the showers 
arising from the decay of such mini black holes and how it compares with a normal (not 
black hole mediated) cosmic ray event, due, for instance, to a primary proton with the 
same energy colliding with an atmospheric nucleon (the "benchmark" event used here). 
The comparison will be based on appropriate observables of the type mentioned above. 

Our incomplete control of the quantum gravity/string theory effects, of the physics 
of low energy non-perturbative QCD and of the nature of the quark-gluon plasma phase 
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in QCD, makes a fully general analysis of the above phenomena impossible at this stage. 
To proceed, we made the following simplifying assumptions and approximations. (1) 
The brane tension was assumed much smaller than the fundamental gravity scale, so it 
does not modify the flat background metric. It is not clear at this point how severe 
this assumption is, since it is related to the "cosmological constant problem" and to 
the concrete realization of the Brane- World scenario. (2) The black hole was assumed 
to evaporate instantly, leading to initial "partons", whose number and distributions are 
obtained semiclassically. No virtual holes were discussed and no back reaction was taken 
into account. (3) The initial decay products were assumed to fly away and hadronize, with 
no intermediate formation of a quark-gluon plasma or of a disoriented chiral condensate 
(DCC). (4) We used standard simulation programs for the investigation of the extensive 
air showers produced in the cases of interest. To this purpose, we have decided to use 
the Monte Carlo program CORSIKA [78] with the hadronic interaction implemented in 
SIBYLL [62] in order to perform this comparison, selecting a benchmark process which 
can be realistically simulated by this Monte Carlo, though other hadronization models are 
also available [83]. Finally, (5) a comment is in order about our selection of benchmark 
process and choice of interesting events. In contrast to the case of a hadronic primary, the 
mini black hole production cross section due to the collision of a > 10 3 TeV neutrino with a 
parton is of the order of the weak interaction neutrino-parton cross section [99]. It would, 
thus, be interesting to compare the atmospheric showers of a normal neutrino-induced 
cosmic ray event to one with a black hole intermediate state. Unfortunately, at present 
neutrinos are not available as primaries in CORSIKA, a fact which sets a limitation on our 
benchmark study. However, it has to be mentioned that neutrino scattering off protons 
is not treated coherently at very high energy, since effects of parton saturation have not 
yet been implemented in the existing codes [27]. As shown in [88] these effects tend to 
lower the cross section in the neutrino case. For a proton-proton impact, the distribution 
of momenta among the partons and the presence of a lower factorization scale should 
render this effect less pronounced. For these reasons we have selected as benchmark 
process a proton-to-air collision at the same depth (X ) and with the same energy as 
the corresponding "signal event". In order to reduce the large statistical fluctuations in 
the formation of the extensive air showers after the collisions, we have chosen at a first 
stage, in the bulk of our work, to simulate collisions taking place in the lower part of the 
atmosphere, up to 1 km above the detector, in order to see whether any deviation from 
a standard scattering scenario can be identified. Another motivation for the analysis of 
such deeply penetrating events is their relevance in the study of the possibility to interpret 
the Centauro events as evaporating mini black holes [99]. A second group of simulations 
have been performed at a higher altitude, for comparison. 
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The present chapter consists of seven sections, of which this Introduction is the first. 
In Section 7.3 we briefly describe the D-brane world scenario, in order to make clear 
the fundamental theoretical assumptions in our study. A brief review of the properties of 
black holes and black hole evaporation is offered here, together with all basic semiclassical 
formulas used in the analysis, with the dependence on the large extra dimensions shown 
explicitly. In Section 7.4 a detailed phenomenological description of the modeling of the 
decay of the black hole is presented, which is complementary to the previous literature 
and provides an independent characterization of the structure of the decay. Incidentally, 
a Monte Carlo code for black hole decay has also been presented recently [89]. We recall 
that this description -as done in all the previous works on the subject- is limited to the 
Schwarzs child phase of the lifetime of the mini black hole. The modeling of the radiation 
emission from the black hole - as obtained in the semiclassical picture - (see [84] for an 
overview) is performed here independently, using semi-analytical methods, and has been 
included in the computer code that we have written and used, and which is interfaced 
with CORSIKA. Recent computations of the greybody factors for bulk/brane emissions 
[84], which match well with the analytical approach of [85] valid in the low energy limit 
of particle emission by the black hole, have also been taken into account. Section 7.5 
contains our modeling of the hadronization process. The hadronization of the partons 
emitted by the black hole is treated analytically in the black hole rest frame, by solving 
the evolution equations for the parton fragmentation functions, making use of a special 
algorithm [25] and of a specific set of initial conditions for these functions [87]. After a 
brief discussion in Section 7.6 of the transformation of the kinematics of the black hole 
decay event from the black hole frame to the laboratory frame, we proceed in Section 7.7 
with a Monte Carlo simulation of the extensive air showers of the particles produced by 
taking these particles as primaries. The simulations are quite intensive and have been 
performed on a small computer cluster. As we have already mentioned, in this work we 
focus on the multiplicities, on the lateral distributions of the events and on the ratio 
of electromagnetic to hadronic energies and multiplicities and scan the entire ultra high 
energy part of the cosmic ray spectrum. Our results are summarized in a series of plots 
and are commented upon in the final discussion in Section 7.8. 

7.3 TeV Scale Gravity, Large Extra Dimensions and 
Mini Black Holes 

The theoretical framework of the present study is the D-brane world scenario [10]. The 
World, in this scenario, is 10 dimensional, but all the Standard Model matter and forces 
are confined on a 4 + ul dimensional hypersurface (the -D3 + „ L -brane). Only gravity with 
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a characteristic scale M* can propagate in the bulk. The til longitudinal dimensions are 
constrained experimentally to be smaller than 0(TeV~ l ). However, for our purposes 
these dimensions may be neglected, since the Kaluza-Klein excitations related to these 
have masses at least of OiTeV' 1 ), too large to affect our discussion below. Consistency 
with the observed Newton's law, on the other hand, leads to the relation Mp h = M™ +2 V n , 
between Mpi ~ 10 19 GeV, the fundamental gravity scale M* and the volume V n of the 
n = 6 — riL dimensional compact or non-compact transverse space. A natural choice for 
M*, dictated a priori by the "gauge hierarchy" puzzle, is M* = 0(M W ) = 0(1 TeV), 
while the simplest choice for the transverse space is an n-dimensional torus with all radii 
equal to R. Thus, one obtains a condition between the number n and the size R{n) of the 
transverse dimensions. Notice that under the above assumptions and for all values of n, R 
is much larger than 10~ 33 cm, the length scale at which one traditionally expects possible 
deviations from the 3-dimensional gravity force, and the corresponding dimensions are 
termed "large extra dimensions" (LED). For n = 2 one obtains R{n = 2) of the order of a 
fraction of a mm. At distances much smaller than R one should observe 3 + n-dimensional 
Newton's law, for instance, as in torsion balance experiments [91]. Current bounds on the 
size of these large extra dimensions and on M* come from various arguments, mostly of 
astrophysical (for instance M* > 1500 TeV for n = 2) or cosmological (M* > 1.5 TeV for 
n — 4) origin [84]. A larger number of LED (n) translates into a reduced lower bound on 
M*. It should be pointed out, that in general it is possible, even if "unnatural", that the 
transverse space has a few dimensions large and the others small. Here we shall assume 
a value of M* of order ITeV, neglect the small extra dimensions and treat the number of 
LED (n) as a free parameter. 

The implications of the existence of LED are quite direct in the case of black hole 
physics. The black hole is effectively 4-dimensional if its horizon (r H ) is larger than the 
size of the extra dimensions. In the opposite case (th <C R, or equivalently for black hole 
masses Mbh <S 10 13 kg for n = 6 [99]) it is 4 + n dimensional, it spreads over the full 
space and its properties are those of a genuine higher dimensional hole. According to some 
estimates, over which however there is no universal consensus [119], black holes should 
be produced copiously [48] [67] in particle collisions, whenever the center of mass energy 
available in the collision is considerably larger than the effective scale M* (y/s » M*). 
With M* ~1 TeV, one may contemplate the possibility of producing black holes with 
masses of order a few TeV. 

Their characteristic temperature Th is inversely proportional to the radius th of the 
horizon, or roughly of order M* and evaporate by emitting particles, whose mass is smaller 
than Th- The radiation emitted depends both on the spin of the emitted particle, on the 
dimension of the ambient space and on the amount of back-scattering outside the horizon, 
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contributions which are commonly included in the so called "greybody factor", which are 
particularly relevant in the characterization of the spectrum at lower and at intermediate 
energy. A main feature of the decaying mini black hole is its large partonic multiplicity, 
with a structure of the event which is approximately spheroidal in the black hole rest 
frame. 

Once produced, these mini black holes evaporate almost instantly. The phenomenolog- 
ical study of 4-dimensional black holes of large mass and, in particular, of their Hawking 
radiation [90] [103] [93], as well as the study of the scattering of states of various spins 
(s = 0,1/2,1) on a black hole background, all performed in the semiclassical approxi- 
mation, have a long history. For rotating black holes one identifies four phases charac- 
terizing its decay, which are (1) the balding phase (during which the hole gets rid of its 
hair); (2) the spin-down phase (during which the hole slows down its rotational motion); 
(3) the Schwarzschild phase (the usual semiclassically approximated evaporation phase) 
and, finally, (4) the Planck phase (the final explosive part of the evaporation process, 
with important quantum gravitational contributions). Undoubtedly, the best understood 
among these phases is the Schwarzschild phase, which is characterized by the emission 
of a (black body) energy spectrum which is approximately thermal, with a superimposed 
energy-dependent modulation, especially at larger values of the energy. The modulation 
is a function of the spin and is calculable analytically only at small energies. Extensions of 
these results to 4+n dimensions are now available, especially in the Schwarzschild phase, 
where no rotation and no charge parameter characterize the background black hole solu- 
tions. Partial results exist for the spin down phase, where the behaviour of the greybody 
factors have been studied (at least for 1 additional extra dimension) both analytically 
and numerically. The Planck phase, not so relevant for a hole of large mass (say of the 
mass of the sun (M ~ 2 x 10 33 gr) which emits in the nano-Kelvin region, is instead very 
relevant for the case of mini-black holes, for which the separation between the mass of 
the hole and the corresponding (effective) Planck mass M* gets drastically reduced as the 
temperature of the hole raises and the back-reaction of the metric has to be taken into 
account. 

In the discussion below we shall use the semiclassical formulas derived for large black 
holes in the Schwarzschild phase and naively extrapolate them to the mini black holes as 
well. This is not, we believe, a severe approximation for the phenomena we shall discuss. 
As the hole evaporates, it looses energy, its mass decreases, its temperature increases and 
the rate of evaporation becomes faster. Thus, the lifetime of the hole is actually shorter 
than the one derived ignoring the back reaction. As we shall see below, the naive lifetime 
is already many orders of magnitude smaller than the hadronization time. This justifies 
the use of the "sudden approximation" we are making of the decay process and explains 
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why the neglect of the back reaction is not severe. 

We recall that the metric of the 4 + n dimensional hole in the Schwarzschild phase is 
given by [101] 



ds 1 



1- — 



r H \n+i 



dt 2 - 



1- — 



r H \n+i 



dr 2 -r 2 d{l 2 2+n , (7.1) 



where n denotes the number of extra spacelike dimensions, and dQ^+n i s the area of the 
(2 + n)-dimensional unit sphere which, using coordinates < ip < 2n and < 0j < 7r, 
with i — 1, n + 1 takes the form 

dn 2 2+n = d9 2 n+1 + sin 2 9 n+1 [d6 2 n + sin 2 9 n ( ... + sin 2 9 2 « + sin 2 9 1 dip 2 ) )) . (7.2) 

The temperature T# of the black hole is related to the size of its horizon by [101] 

T„ = £±! (7.3) 
47T r# 

and the formula for the horizon r# can be expressed in general in terms of the mass of 
the black hole Mbh and the gravity scale M* [101] 

r „ = w^-(!£M^ (74) 

v 7 ^* V M, y ^ n + 2 J { ' 

For n = and M* = M Pi ~ 10 19 GeV it reproduces the usual formula for the horizon 
(r H = 2GM BH ) of a 4 dimensional black hole. For n > the relation between r H and 
Mbh becomes nonlinear and the presence of M* in the denominator of Eq. (7.4) in place 
of Mpi increases the horizon size for a given M BH . For M B ^/M* ~ 5 and M* = 1 TeV 
the size of the horizon is around 10~ 4 fm and decreases with increasing n. 

In the Schwarzschild/spin-down phase, the number of particles emitted per unit time 
by the black hole as a function of energy is expressed in terms of the absorption/emission 
cross sections a^ n {uo) (or equivalently of the greybody factors T(uj)), which, apart from n, 
depend on the spin (s) of the emitted particle, the angular momentum (j) of the partial 
wave and the corresponding energy (cu), 

dNU(u) x . jgH oo 2 



dtdw y 2ir 2 exp (u/T H ) ± 1 



Multiplying the rate of emitted particles per energy interval dN^ s \u) / dtdu by the particle 
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energy u one obtains for the power emission density 

dtdu y 2vr 2 exp(u/T H )±l U [ ' 

where the sum is over all Standard Model particles and the +(— ) in the denominator 
correspond to fermions (bosons), respectively. a^ n are the cross sections for the various 
partial waves and depend on the spin s of each particle. We recall, that in the geometric 
optics approximation a black hole acts as a perfect absorber of slightly larger radius r c 
than th [HO], which can be identified as the critical radius for null geodesies 



- n + 3\V(n+i) i n + 3 



The optical cross section is then defined in function of r c (or equivalently Th via 
Eq. (7.7)), such that A k , the effective surface area of the black hole hole projected over a 
k-dimensional sub-manifold becomes [56]: 

2 . . fc-2 



d- 1 \ d - [i d- 1 



2 



1 !-7-! [^ij r k H 2 (7.8) 



and 



2n 2 



is the volume of a k-sphere. 



It is convenient to rewrite the greybody factors as a dimensionless constant T s = a s j ' A± 
normalized to the effective area of the horizon At, obtained from (7.8) setting k — 4 and 
d = 4 + n 

'n + 3\ 2/(n+1) n + 3 2 



Ai = 4 * {—) —i r «' (7 - 10) 

and replacing the particle cross section a in terms of a thermal averaged graybody factor 
ri(ri/2 = 2/3, Ti = 1/4, r = 1, % denoting the spin or species [102]). Eqs. (7.5) integrated 
over the frequency give (for particle %) 



^ = a(n,r H )n (7.11) 

with 

a(n, r H ) = A T t T(3) C(3) c t A 4 T 3 H , (7.12) 
where q is the number of degrees of freedom of particle i and fi is defined by the integral 
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(sj is the spin) 



I rf V/^-(-i^ =/ - r(3)c(3) ^ (7 - 13) 



from which /j = 1 (/, = 3/4) for bosons (fermions). These numbers depend on the 
dimension of the brane, which in our case is 3. T(x) and ((x) are the Gamma and the 
Riemann function respectively. Since A 4 depends on the temperature (via r#), after some 
manipulations one obtains 

1 /„ + <*\ 2/(n+l) 

A 4 T| = — (n + 3)(n + l)T fl (7.14) 



4tt V 2 

and 



§ = 8^ 2V0. + .» (» + Dr(3)C(3)r,cr ff . (7.15) 
Summing over all the particles i we obtain the compact expression 

dN ^(j2fiT iC ^r(3)((3)T H (7.16) 



dt 2vr 
with 

_ r t (n+l)(n + 3)(" +3 )/(" +1 ) 
' ~~ 47T 2 2 2 /(™ +1 ) 

The emission rates are given by 



(7.17) 



( , 18) 



J», « 4 x 3.7 x 10- i!i^__-XiZ ^ r . , (7 . 19) 



(yl + 3)(n+3)/(n+l)( n+ l) / T H 

2 2 /("+i) ^GeV 

^ 4xia>lf M^!M (|l) .-, ( , 20) 

for particles with s = 0, 1/2, 1, respectively. Furthermore, integration of Eq. (7.6) gives 
for the black hole mass evolution 

dM dE . 

~dT = -^ = ^ Th)Th 

= ^(Z/^c,) T(4)C(4)T|, 

(7.21) 



with 



= 2^Ete r */'i) ^r(4)c(4) 



(7.22) 
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where now f[ = 1 (7/8) for bosons (fermions). Taking the ratio of the two equations 
(7.21) and (7.16) we obtain 

dN M 1 



dM \(3) T H 
4ir6(n) 

n+ 1 M* \M, 



= P^TTT TT) . (7-23) 



where we have defined 



and 



'< n >=i-^rJ 7i' (7 - 24) 

n E,Q/,r t r(3)C(3) 

Ei Q /'i I\ r(4) c(4) ' (7 - 25) 

This formula does not include corrections from emission in the bulk. 

In the "sudden approximation" in which the black hole decays at its original formation 
temperature one easily finds N = ^^p-^, where E is the energy spectrum of the decay 
products, and using Boltzmann statistics iV £s one obtains the expression [48] 

This formula is approximate as are all the formulas for the multiplicities. A more accurate 
expression is obtained integrating Eq. (7.23) to obtain 

and noticing that the entropy of a black hole is given semiclassically by the expression 



So 



n + lM 
n + 2T 

n+2 



n + 2 V M, 



0(n), 



one finds that 



N = P S , (7.28) 



which can be computed numerically for a varying n. As one can see from Fig. 7.1, the 
two formulas for the multiplicities are quite close, as expected, but Eq. (7.26) gives larger 
values for the multiplicities compared to (7.28) as noted by [33]. Other expressions for 
the multiplicities can be found in [34]. 
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Figure 7.1: Multiplicities computed with Eq. (7.28) (A) and Eq. (7.26) (B) for a varying 
number of extra dimensions n. 



Since the number of elementary states becomes quite large as we raise the black hole 
mass compared to the (fixed) gravity scale, and given the (large) statistical fluctuations 
induced by the formation of the air shower, which reduce the dependence on the multi- 
plicity formula used, we will adopt Eq. (7.26) in our simulations. Overall, in the massless 
approximation, the emission of the various species for a 3-brane is characterized by ap- 
proximately 2% into spin zero, 85% into spin half and 13% into spin one particles, with 
similar contributions also for the power emissivities. These numbers change as we vary the 
dimension of the brane (d) and so does the formula for the emissivities, since the number 
of brane degrees of freedom (q(c?)) has to be recomputed, together with the integrals on 
the emission spectra (fi(d)) [1]. 

The integration of the equation for the power spectrum, in the massless approximation, 
can be used to compute the total time of decay (assuming no mass evolution during the 

decay) 

1 /M BH \ (n+3)/(n+1) 



r ~sihrJ < 7 - 29 ' 

which implies that at an energy of approximately 1 TeV the decay time is of the order 
of 10~ 27 seconds. Therefore strong interaction effects and gravity effects appear to be 
widely separated and hadronization of the partons takes place after their crossing of 
the horizon. The black hole is assumed to decay isotropically (s-wave) to a set of N 
elementary states, selected with equal probability from all the possible states available in 
the Standard Model. We mention that in most of the analysis presented so far [60, 7, 1] 
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the (semiclassical) energy loss due to bulk emission has not been thoroughly analyzed. 
We will therefore correct our numerical studies by keeping into account some estimates 
of the bulk emission. 

7.4 Modeling of the Black Hole Decay 

The amount of radiation emitted by the black hole in the ED is viewed, by an observer 
living on the brane, as missing energy compared to the energy available at the time when 
the black hole forms. From the point of view of cosmic ray physics missing energy channels 
imply reduced multiplicities in the final air shower and modified lateral distributions, these 
two features being among the main observables of the cosmic ray event. However, since 
the initial energy of the original cosmic ray is reconstructed by a measurement of the 
multiplicities, an event of reduced multiplicity will simply be recorded as an event of 
lower energy. It is then obvious that an additional and independent reconstruction of the 
energy of the primary cosmic ray is needed in order to correctly identify the energy of 
these events. 

In our study we will compute all the observables of the induced air shower using 
both the lab frame (LF) and the black hole frame (BHF) to describe the impact and the 
formation of the intermediate black hole resonance. Also, in the simulations that we will 
perform, the observation level at which we measure the properties of the air showers will be 
selected to take properly into account the actual position of a hypothetical experimental 
detector. The target of the first impact of mass M is assumed to be a nucleon (or a 
quark) at rest in the atmosphere and the center of mass energy, corrected by emission 
loss in the bulk, is made promptly available for an instantaneous black hole formation 
and decay. We will also assume that the energy E\ of the incoming primary varies over 
all the highest part of the cosmic ray spectrum, from 10 15 eV up to 10 20 eV. 

We denote by (3 the speed of the black hole in the lab frame. In our notations, E* is 
the typical energy of each elementary state in the decay products (parton, lepton) in the 
BH frame and P* is its corresponding momentum. 

We will assume that a black hole decays "democratically" into all the possible partonic 
states, proportionally to the number of Standard Model states which are available to it 
at a given energy. 

The energy per partonic channel will be appropriately weighted and we will assume 
that each parton (/) will decay into a final state hadron h (carrying a fraction x of the 
original momentum) , with a probability distribution given by the corresponding fragmen- 
tation function Df(x,Q 2 ), which is evolved from some low energy input scale Qo up to 
the relevant scale characterizing the decay. This is given by the available energy per 



196 



7.4. Modeling of the Black Hole Decay 



fundamental state, equally distributed among all the states. 

The quantification of the injection spectrum involves a computation of the relevant 
probabilities for the formation of all the possible hadronic/leptonic states prior to the 
simulation of the air shower. Let's briefly elaborate on this. 

To move from the parton level to hadron level, we let D^(x,Q 2 ), D^(x,Q 2 ), and 
Dg(x, Q 2 ) be the fragmentation functions of N F quarks q, antiquarks q, and of the gluon 
g, respectively, into some hadron h with momentum fraction x at the scale Q. From 
the fragmentation functions we obtain, for each hadron h, the mean multiplicity of the 
corresponding s-wave and the corresponding average energy and momentum. Specifically 
we obtain 

<D h > = W* dzD h f (z,Q 2 ) (7.30) 

for the probability of producing a hadron h, and 

E* h = W 1 zdzD h f (z,Q 2 ) (7.31) 

for the average energy of the same hadron. We recall that z m i n is the minimal fraction 
of energy a hadron (h), of mass m h , can take at a scale Q, and can be defined as z min = 
rrif l /(Q/2). In practical applications one can take the nominal value z m i n = 0.05 for every 
hadron, without affecting much the mean multiplicities and the related probabilities. This 
implies that 

< D h r > +]T < D) > + < D h g > + < >= Pv h (7.32) 

/ 

together with the condition J2h P r /i — 1) where the sum runs over all the types of hadrons 
allowed by the fragmentation. In all the equations above, the fragmentation takes place 
at the typical scale Q = E/N, scale at which the moments are computed numerically. 
For the identification of the probabilities it is convenient to organize the 123 fundamental 
states of the Standard Model into a set of flavour states ((?/), with / running over all the 
flavours except for the top quark, where in (gy) we lump antiquark states and color states, 
plus some additional states. The weight of the (<?/) set is pf = 2 x 2 x 3/123, where the 
factors 2 and 3 refer to spin, quark-antiquark degeneracy and color. It is worth to recall 
that quark and antiquark states of the same flavour have equal fragmentation functions in 
all the hadrons, and this justifies the q/q degeneracy of the set. The additional states are 
the gluon (g) with a weight p g = 2 x 8/123, the photon (7), with a weight p 7 = 2/123 and 
the remaining states (r) in which we lump all the states which have been unaccounted for, 
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whose probabilities p r are computed by difference. These include the top and the antitop 
(12/123), the W's and Z (9/123) and the leptons (24/123). The fragmentation functions 
into hadrons, corresponding to these states, < > are computed by difference from the 
remaining ones < D h g >, < Dj > and < >, which are known at any scale Q from 
the literature. Beside the favour index f — u,d, c, s, b, introduced above, we introduce a 
second index i running over the (r) states, the photon and the gluons (i = g,j,r). 



The probability of generating a specific sequence of N states in the course of the 
evaporation of the black hole is then given by a multinomial distribution of the form 

f(n f ,n t ,p f , Pi ) = N ; , 11^11^ (7-33) 
11/ H f- ili <H- f i 

which describes a typical multi-poissonian process with N trials. Notice that, to ensure 
proper normalization, we need to require that 

[ J //,! = n 9 !n 7 !n r ! 

% 

= ngln^N-^nf-ng-n^l (7.34) 
/ 



The computation of the cumulative probabilities to produce any number of hadrons of 
type h by the decay of the black hole are obtained from the multinomial distribution mul- 
tiplied by the fragmentation probabilities of each elementary state into (h) and summing 
over all the possible sequences 

Prcum h - E n n^nn' U ^ < ^ ^ U ^ < ^ ^ ' ^ 

n f ,ni 11/ u f- Hi 'H- / j 

A possible way to compute Pr cum ^ when is large is to multiply the multinomial 
distribution by a suppression factor Exp[— A(J2i n i + J2f n f ~ N)], with A a very large 
number, and interpret this factor as a Boltzmann weight, as in standard Monte Carlo 
computations of the partition function for a statistical system. Simulations can be easily 
done by a Metropolis algorithm and the configurations of integers selected are those for 
which the normalization condition A^ = J2i n i + J2f n f is satisfied. In our case, since we are 
interested only in the mean number of hadrons produced in the decay and in their thermal 
spectrum, the computation simplifies if we average over all the relevant configurations. 
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7.5 Fragmentation and the Photon Component 



The evolution with Q 2 of the fragmentation functions is conveniently formulated in terms 
of the linear combinations 



N F 



1=1 



D^(x,Q 2 ) = D%.(x,Q 2 ) + D^{x,Q 2 ) - —D*(x,Q 2 ) , 



Dl lt (x,Q 2 ) = D h qi (x,Q 2 )-Dl.(x,Q 2 ), 



(7.36) 
(7.37) 
(7.38) 



as for these the gluon decouples from the non-singlet (+) and the asymmetric (-) combi- 
nations, leaving only the singlet and the gluon fragmentation functions coupled; 



d 



Q2 d& Dh ^ Q2) = I P(+) {*s{Q 2 )) ® D h {+) ,{Q 2 )] (x) , 



Q2 dQi D ^ {x > Q2) = [ P (-)(*sm)®D^ t (Q 2 )](x), 



Q 2 -^D*(x,Q 2 ) = [P*(a s {Q 2 ))®D*{Q 2 



x 



+2N, 



P G ^ q (a s (Q 2 )) ® D*(Q 2 )] (x) 
P 9 ^ 9 {* S {Q 2 )) D h g (Q 2 )} (x) . 



+ 



The kernels that appear in the equations above are defined by 



(7.39) 
(7.40) 



P q ^ G (a s (Q 2 )) ® D h G (Q 2 )} (x) , (7.41) 



(7.42) 



P (+) (x, a s (Q 2 )) = P^ q (x, a s (Q 2 )) + P g % (x, a s (Q 2 )) , (7.43) 
Pe (x, a s (Q 2 )) = P (+) (x, a s (Q 2 )) + 2N F P q t q (x, a s (Q 2 )) , (7.44) 



P ( _) (x, a s (Q 2 )) = P^ q (x, a s (Q 2 )) - P g % (x, a s (Q 2 )) , (7.45) 
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with a s (Q 2 ) being the QCD coupling constant. In the perturbative expansion of the 
splitting functions, 



P(x,a s (Q 2 )) = ^(x) + [^) 2 P {1 \*) + O 



'0is{Q 

2ir 



2\\ 3 N 



The timelike kernels that we use are given by 

"3 



P^(* 



p(0) / 

q^q V 



x 



C f 



-5(1 -x) + 2 







(l-x) 


— 1 — X 




f 



= ^(x) = 0, 

1 + (1 -x) 2 



2NpTjt x 2 + (1 - 



+2iV c 



1 



1 — x/ + x 
The formal solution of the equations is given by 



H 2 + x-x 2 



(7.46) 



(7.47) 
(7.48) 
(7.49) 
(7.50) 

(7.51) 



rlog{Q 2 /Ql) n (r>2\ 

D h a (*,Q 2 ) = D h a (x,Ql) + J o d\og{Q*/Ql)^^Y<[P^MQ 2 ))®D^ 

(7.52) 

where Qo is the starting scale of the initial conditions, given by D(x,Qq). At leading 
order in a s , we solve this equation using a special ansatz 



Z^,Q 2 )=£^l g^ (g2) ^ 



a,(Qo) 



(7.53) 



and generating recurrence relations at the n+l—th order for the A n+1 coefficients in terms 
of the A n [25]. It is easy to see that this corresponds to the numerical implementation of 
the formal solution 

D h f (x, Q 2 ) = Exp (tP®) D h f (x, Ql) (7.54) 

with t = (a s (Q 2 )/a s (Ql)), where the exponential is a formal expression for an infinite 
iteration of convolution products. We show in Figures 7.2-7.4 results for some of the 
fragmentation functions into pions, kaons and protons computed for a typical parton 
scale of 200 GeV. 

The photon contributions to the decay of the black hole is treated separately. The 
evolution equation for the fragmentation functions of photons and parton fragmentation 
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Figure 7.2: Fragmentation functions into 7r =l= at 200 GeV 
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Figure 7.3: Fragmentation functions into K ± at 200 GeV. 
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Figure 7.4: Fragmentation functions into p/p at 200 GeV. 
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into photon D^(x, Q 2 ), D^(x, Q 2 ) satisfy at leading order in a ern (the QED fine structure 
constant) and a s (the QCD coupling), the evolution equations [86] 

y '^ 1 '*^'^ <7 - 55) 

which can be integrated with the initial conditions D^(x, Q 2 ) = 5(1 — x), and 



dDVx,Q' 



a 



P^(x)(g)DVx,Q 2 ) (7.56) 



dlnQ 2 2tt 
which can also be integrated with the result [65] 

D,/ g (x, Q 2 ) = ^W*) In |j + Ql). (7.57) 

In Eq. (7.57) the second term is termed the hadronic boundary conditions, which come 
from an experimental fit, while the first term is the pointlike contributions, which can be 
obtained perturbatively. In [3] a leading order fit was given 

ry r- 

DUx, Ql) = - P q ^(x) ln(l - x) 2 - 13.26 (7.58) 

at the starting scale Q = 0.14 GeV. The kernels in these cases are given by simple 
modifications of the ordinary QCD kernels, for instance P g _> 7 (x) = e 2 /CFPq-> g {x) [86]. 

For the fragmentation function of quarks to photons with virtuality M 7 , the pertur- 
bative result is given in [105] 



D, /q (x,Q 2 ,M 2 ) = e 2 



a 
2^ 



l + (l-x) 2 . xQ 2 / P 2 
In — — - — x 1 — 



x \ xQ 2 



(7.59) 



where P 2 is the virtuality of the photon. The gluon to photon transitions are neglected, 
since P g ^ vanishes in leading order. 

We recall that each elementary state emitted is characterized by an average energy 
given by (e) = M BH /(N). 

The leptonic component e ± , /x 1 * 1 , produced by the decay is left unaltered and provides 
an input for the air shower simulator as soon as these particles cross the horizon. The r ± s 
are left to decay into their main channels, while the hadronization of the u, d, s, c quarks 
and the gluons is treated with our code, that evolves the fragmentation functions to the 
energy scale (e). The top (t) quark is treated consistently with all its fundamental decays 
included; hadronization of the b quark is treated with suitable fragmentation functions 
and also involves a suitable evolution. As we vary Mbh and we scan over the spectrum of 
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Table 7.1: Initial conditions. For each couple of parton and hadron, the upper num- 
ber in the box is the probability for the parton / to hadronize into the hadron h, 

(XL„ D f( z ' Q"> dz ) I (-tin ' while the lower number is the average energy 

fraction of h, J^ m _ n zD l j(z,Q)dz. In this table the energy of the parton is Q = 1.414 GeV 
for u, d, s andfl^Q = 2m c = 2.9968 GeV for c and Q = 2m b = 9.46036 GeV for the b 
quark, generated via the set of ref. [87] 

the incoming cosmic rays the procedure is repeated and rendered automatic by combining 
in a single algorithm all the intermediate steps. Tables 1 and 2 contain the results of a 
renormalization group analysis of the fragmentation functions for all the partons (except 
the top quark), where we show both the initial conditions at the input scale, whose lowest 
value is Q — 1.414 GeV, and the results of the evolution, at a final scale of Q — 200 GeV, 
the initial set being taken from ref. [87]. 

This concludes the computation of the probabilities for each hadron/lepton present 
in the decay products of the mini black hole. It is reasonable to assume that these 
particles will be produced spherically, since higher angular momenta are suppressed by 
the corresponding centrifugal barrier. However, the analysis of the shower profile has to be 
performed in the lab frame. This requires the transformation of the initial configurations 
above to the laboratory frame, which is exactly what is discussed next. 

7.6 Sphericity and Boost 

The transformation from the black hole frame (BHF) to the laboratory frame (LF) is 
performed by a Lorentz boost with speed (3, the speed of the black hole in the LF. 
Assuming that the black hole is produced in the collision of a primary of energy E 1 in 
the LF and negligible mass compared to E±, with a parton of mass M in the atmosphere, 
one obtains (3 = E\j{E\ + M). A spherical distribution of a particular particle of mass m 
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Table 7.2: For each couple of parton/hadron, the first number is the probability of frag- 
mentation of the parton / into the hadron (jl min D l j(z,Q)dz^j / J2h' {jl min Df'(z, Q)dz^j , 
while the second is the average energy fraction of h, l] mm zD h f (z, Q)dz. The energy of the 
parton is Q = 200 GeV. 



among the decay products in the BHF is transformed to an elliptical one, whose detailed 
form is conveniently parametrized by 



P 1 



M 



g" = ^= Sl±i< (7.00) 

V (l-fl)" 2 

where (3* = P*/E* is the speed of this particle in the BHF, the ratio of its BHF momentum 
and energy. Figure 7.5 depicts the relevant kinematics. A particle emitted in the direction 
9* in the BHF, is seen in the direction 9 in the LF, with 

1 sin 9* 

tanfl = 0^— -. (7.61) 

For values of g* > 1 the shape of the 1-particle distribution in the LF is characterized by 
a maximum angle of emission 



t&n9 max \ = (7.62) 



which for g* — 1 is equal to 90°. Only for g* < 1 there is backward emission in the LF. As 
a relevant example, let us consider the case of a hadron of mass m = 1 GeV and energy 
E* = 100 GeV emitted by a black hole, formed by an initial primary of energy E 1 = 1000 
TeV, which hit a quark of mass M ~ lOMeV to form a black hole. The gives g* = 1.00005 
and the corresponding maximum angle in the LF is 9 max — t&ri9 max = 1.4 x 10~ 2 , giving 



204 



7.6. Sphericity and Boost 




an angular opening of the decay products of about 2 degrees. 

As shown in Figure 7.5, for g* > 1, which is relevant for our purposes, the mapping 
from 9* to 9 is not one to one. In a given direction 9 in the LF, one receives particles 
emitted in two different directions 9*^ in the BHF. They satisfy [24] 



dcos9 ' \P*J (±cos6VK) 

where 



7(1 — v 2 cos 2 9) 

and 



1/2 



(7.63) 



K= l + 7 2 (l-#* 2 )tan 2 fl, (7.64) 
and with the momenta P ± and energies E ± of the two branches given by 

. P*cos9(g*±VK) 



(7.65) 



m ( 7* ± v cos 9 (v* 2 Y 2 - v 2 ^ 2 sin 2 9) 1 

E± = — 7T 2 2 m — ( 7 - 66 ) 

7(1 — v z cos z 9) 

respectively. In the above formulas 7 _1 = ^/l — (3 2 , v(v*) is the speed of the hadron in 
the LF (BHF), and (7*) _1 = Tn/E^ = y/l — v* 2 . For massless final state particles, in 
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particular, these relations become 

P = E= P * (7.67) 
7(1 — cos 9) 

and reduce to the familiar Doppler formula when 9 = 0. 

The probability distribution Wj*(cos0*, 0) of a hadron (h) as a function of the direction 
Q = (cos#*,0) in the BHF, assumed spherically symmetric and normalized to the total 
probability Pr h of detecting this hadron among the decay products with TV elementary 
states, is 

W* h (cos9*) = ^. (7.68) 
The corresponding one in the LF is 

W h (cos9) = 2^Vw fc *(cose* ± ). (7.69) 
_i_ a cos (7 

In the special case g* — 1, the probability distribution simplifies to 

cos 

W h {cos 9) = 2Pv h — -— , (7.70) 

7 2 (1 — p 2 cos 9 2 ) 2 

peaked in the forward direction, symmetric around the maximum value, obtained for 
9 = and equal to 27 2 , while the momentum distribution is 

7(1 — p l cos9 2 ) 

As we have already mentioned, the structure of the partonic event (and, similarly, of 
the hadronic event after fragmentation) is characterized by the formation of an elliptical 
distribution of partons, strongly boosted toward the detector along the vertical direction. 
Each uniform (s-wave) distribution is strongly elongated along the arrival direction (due 
to the large speed of the black hole along this direction) and is characterized by two 
sub-components (W^), identified by a ± superscript. Their sum is the total probability 
distribution given in (7.63). The "+" momentum component is largely dominant and 
strongly peaked around the vertical direction with rather small opening angles and this 
behaviour can be analyzed numerically with its iV dependence. In the explicit identifica- 
tion of the two independent distributions ± in terms of the opening angle 9, as measured 
in the LF, we use the relations 

W ± {9) = -W ± (cos9*) | dcos ^ 9 *J ± | sin g ? (7.72) 
w 2 v 71 dcos9 1 y ' 
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where we have introduced a factor of 1/2 for a correct normalization of the new distribution 
in the 9 variable. In Figures 7.6 and 7.7 we show the structure of these distributions in 
the LF. Both are characterized by a very small opening angle (9) with respect to the 
azimuthal direction of the incoming cosmic ray, W + being the dominant one. Two similar 
plots (Figures 7.8 and 7.9) illustrate the two components P ± as functions of the same 
angle. 

One can easily check that we can now integrate symmetrically on both distributions 
to obtain the correct normalization (to Pr/J for a given hadron (or parton) 

/8max , rQi* r—0i* 

W ± (9)d9 = / ^ W + (9)d9 + / ^ W-{9)d9 + / W~(9)d9 = Pv h (7.73) 

or, equivalently, using cos 9 as a distribution variable 

d ens 9* ± 

W ± (cos9) = W ± (cos9*) | | sintf (7.74) 

d cos 9 

with 

f 1 W ± (cos 9)d cos 9 = f 1 W + (cos 9*)d cos9* + T W- (cos 9*)d cos 9* = Pr h 

J CO S0 'max J COS 01* J — 1 

(7.75) 

and 9*i obtained from (7.62). 
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Figure 7.6: Plot of the W + branch of the probability distribution for an incoming energy 
of the cosmic ray E\ = 10 8 GeV and for various values of N of the elementary partonic 
states emitted during the decay of the black hole 



0.0018 
0.0016 
0.0014 
0.0012 
0.001 
0.0008 
0.0006 
0.0004 
0.0002 




-15 



-10 



E 1= 10® GeV,.N=10 
E 1= 10, GeV,! 
E 1= 10 8 GeV, 



iN=12 
|N=15 



10 



15 



Figure 7.7: Same plot as above but for the W branch of the distribution 
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Figure 7.8: Plot of the P + branch of the momentum distribution in eV for an incoming 
energy of the cosmic ray E\ = 10 8 GeV and for various values of N of the elementary 
partonic states emitted during the decay of the black hole 





Figure 7.9: Same plot as above but for the P branch of the momentum distribution 



Chapter 7. Searching for extra dimensions 



209 



n 


L 

J n 


n 

U 


D 7071 1 


i 
i 


U.UUJOO 


9 




3 


0.62057 


4 


0.60696 



Table 7.3: Fraction f n of i?cM that is bound into the black hole as a function of the 
number n of extra-dimension in head-on collisions. 

7.7 Air Shower Simulations 

The simulation of the events is performed at the last stage, using an air shower simulator. 
We have used CORSIKA [78] with appropriate initial conditions on the spectrum of the 
incoming particles in order to generate the full event measured at detector level. In most 
of the simulations we have assumed that the first impact takes place in the lower part 
of the atmosphere, not far from the level of the detector, at a varying altitude. The 
reason is that one of our interests is the investigation of the possibility that the Centauro 
events may be related to evaporating mini black holes, formed by the collision of weakly 
interacting particles (e.g. neutrinos) which penetrate the atmosphere. Of course, we have 
simulated events happening at higher altitudes as well. We have performed two separate 
sets of simulations, the first set being benchmark events with an "equivalent" proton 
replacing the neutrino-nucleon event, colliding at the same height, the second being the 
signal event, i.e. the black hole resonance. The difference between the first and the second 
set is attributed to the different components of the final state prior to the development 
of the air shower. 

We compute the average number of particles produced in the process of BH evaporation 
using the formula 

M BH = E CM fn (7.76) 

where Ecu is the energy in the center of mass frame in the neutrino-nucleon collision and 
f n is the fraction of E C m that is bound into the black hole as a function of the number 
n of extra-dimensions. Numerical values for /„ in head-on collisions are taken from Ref. 
[50] and reported in Table 7.3. 

The overall shower, defined as the superposition of the various sub-components, de- 
velops according to an obvious cylindrical symmetry around the vertical z-axis near the 
center. We assume in all the studies that the incoming primary undergoes a collision with 
a nucleon in the atoms of the atmosphere at zero zenith with respect to the plane of the 
detector. 

The model of the atmosphere that we have adopted consists of N 2 , 2 and Ar with the 
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Table 7.4: Parameters of the U.S. standard atmosphere. 
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Figure 7.10: Pressure versus altitude for the U.S. standard atmosphere model. 



volume fractions of 78.1%, 21.0% and 0.9% [122]. The density variation of the atmosphere 
with altitude is modeled by 5 layers. The pressure p as a function of the altitude h is 
given by 

p(h) = aj + biexp(-h/ci), i = l,...,4 (7.77) 
in the lower four layers and by 

p(h) =a 5 -h— (7.78) 

c 5 

in the fifth layer. 

The Oj, bi, Ci parameters, that we report in Table 7.4, are those of the U.S. standard 
atmospheric model [78]. The boundary of the atmosphere in this model is defined at 
the height 112.8 km, where the pressure vanishes. In Figure 7.10 we show a plot of the 
pressure {p — X v , also called vertical depth) as a function of the height. 
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This is defined via the integral 



roo 

X v = / p(h')dh' 

J h 



(7.79) 



of the atmospheric density p{h) for zero zenith angle, while the corresponding slant depth 
is given by 

roo ( 1 /2 \ 

(7.80) 



roo / 1/2 N 

X = / pllcosO + -—sm 2 9 
Ji \ 2 Rt j 



for a zenith angle 6 and Rt is the radius of the earth. 

To put into perspective our Monte Carlo study it is convenient to briefly summarize 
the basic features of the theory of cascades on an analytical ground. The theory consists 
of the system of transport equations [64, 49] for the numbers N n (E,X) of particles of 
type n with energy E at height X 



dN n (E,X) N n (E,X) 
dX X n (E) 



CT n 7p A ir 



-N n {E,X) , 



(7.81) 



where X n (E) is their interaction length, r n is their lifetime and 7 the Lorentz-factor 
corresponding to their given energy. In the simple case of an isothermal atmosphere 
PAir = Po ex P( — h/h ) = X/h , at a scale height h 



dN n (E,X) N n (E,X) 1 



dX \ n {E) 
where d n is their decay length, defined by 



d 7 



N n (E,X) 



(7.82) 



d r . 



mc 2 ho 
Ect„X 



(7.83) 



Particles produced at higher energies are also accounted for by an additional term in the 
cascade 



dN n (E,X) 
dX 



1 1 

+ 



-N n (E,X 

+ 2/ N m (E',X) 



X n (E) d n (E)_ 

W mn (E',E) 



(7.84) 



\ m (E>) 
dE' , 



+ 1J¥) D ^ E '- E \ 

describing the change in the number of particles of type n due to particles of type m 
by interaction or decay, integrated over an allowed interval of energy. The functions 
W mn (E', E) are the energy-spectra of secondary particles of type n in a collision of particle 
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m with an air-molecule, while D mn (E',E) are the corresponding decay- functions. The 
advantage of a transport equation compared to a Monte Carlo is, that it provides a rather 
simple analytical view of the development of the cascade across the atmosphere. Most 
common in the study of these equations is to use a factorized ansatz for the solution 
N(E,X) = A(E)B(X), which assumes a scaling in energy of the transition functions 
[64]. In our case an analytical treatment of the cascade corresponds to the boundary 
condition 



with Pr n (E) being the probability that the black hole decays into a specific state n. As we 
have already discussed above, these decays are uncorrelated and Eq. (7.85) is replicated 
for all the elementary states after hadronization. The emission probabilities Pv n (E) have 
been computed by us for a varying initial energy E = fMsu using renormalization 
group equations as described before, having corrected for energy loss in the bulk. The 
interactions in the injection spectrum of the original primaries at our X (X = 517 
g/cm 2 ) has been neglected since this is not implemented in CORSIKA. The showers have 
been performed independently and the results of the simulations have been statistically 
superimposed at the end with multiplicities computed at detector level (X\ = 553 g/cm 2 ). 
We have kept the gravity scale M* constant at 1 TeV and varied the mass of the black 
hole according to the available center of mass energy E. As we have already discussed 
in the previous sections, as a benchmark process we have selected a proton-air impact at 
the same X with the boundary condition 



which occurs with probability 1 (Pr p = 1). 

We are interested both in the behaviour of the multiplicities and in the lateral distri- 
butions of the cascades developed at detector level. For this purpose we have defined the 
opening of the conical shower after integration over the azimuthal angle, as in [27], and 
given the symmetry of the event, we plot only the distance from the center as a relevant 
parameter of the conical shower. 

A varying number of extra dimensions n — 0, 1, . . . , 4 implies a different ratio for bulk- 
to-brane energy emission, a different average number of elementary decaying states and 
different energy distributions among these. We have varied the energy E\ of the primary, 
thereby varying the mass of the black hole resonance, in the interval 10 15 — 10 20 eV. The 
hadronization part of our code has been done by changing the Pi h(E) obtained from a 
numerical solution of the fragmentation functions separately for each value of the energy 
shared. 



N n (E,X ) = Pv n (E) 5{E - fM BH /(N)) 



(7.85) 



N p (E,X ) = Pr p 6(E-fM BH ) 



(7.86) 
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Figure 7.11: Benchmark event: Multiplicities of photons, e ± , /r* 1 at an observation level 
of 4500 m as a function of the altitude of the first impact. 

• Preliminary studies 

We start with the numerical study of the partial and total multiplicities of the 
various sub-components and of their lateral distributions in the benchmark event. 
The results for photons and leptons are shown in Figures 6.11 to 7.14 as functions 
of the altitude of the impact and for two different observation points at 4,500 m 
and 5,000 m, approximately the altitudes of the detectors at Pamir and Chacaltaya, 
respectively. These preliminary plots, based on actual simulations of a proton- 
to-air nucleus impact at 10 15 eV, show a steady growth of the multiplicities of 
the secondaries as we raise the point of first impact above the detector. The 
statistical errors in the simulations (80 uncorrelated events have been collected per 
point) are rather small, quite uniformly over all the altitudes of the impact, and 
indicate a satisfactory stability of the result. The positions of the detectors do not 
seem to have an appreciable impact on the characteristics of the secondaries. As for 
the lateral distributions we observe an increase in the opening of the showers with 
the event altitude, which is more enhanced for the muonic component and for the 
photons and less for the electrons and positrons. Also in this case the statistical 
fluctuations are rather small. On the basis of these results we have selected for the 
remaining simulations a first impact at 5,500 m and the observation (detector) level 
at 5,000 m. However, for comparison, we will also show later the results of a second 
set of simulations that have been performed with the first impact at 15,000 m. 

• Choice of scales and corrections 
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Figure 7.12: As in Figure 7.11, but at an observation level of 5000 m. 
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Figure 7.13: Average radial opening of the shower of photons, e ± , /r* 1 at an observation 
level of 4500 m as a function of the altitude of proton's first interaction. 
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Figure 7.14: As in Fig. 7.13, but at an observation level of 5000m. 

To compare standard and black holes events, we have selected a gravity scale of 
M* = ITeV and varied the black hole mass, here taken to be equal to the available 
center of mass energy during the collision. Therefore, a varying E-y is directly related 
to a varying Mbh and we have corrected, as explained above, for the energy loss 
into gravitational emission. Unfortunately, this can be estimated only heuristically, 
with bounds largely dependent on the impact parameter of the primary collision. 
A reasonable estimate may be of the order of 10 — 15% [85]. Corrections related 
to emission in the bulk have also been included, in the way discussed in previous 
sections. 

• Energy ratios: electromagnetic versus hadronic 

Not all observables are statistically insensitive to the natural fluctuations of the air 
showers. In the study of black hole versus standard (benchmark) events, the study 
of the ratios -/V em /iV hadron and E em / E-^^ron have been proposed as a way to distin- 
guish between ordinary showers and other extra-ordinary ones. Centauro events, 
for instance, have been claimed to be characterized by a rather small ratio of elec- 
tromagnetic over hadronic energy deposited in the detectors, contrary to normal 
showers, in which this ratio is believed to be -E'em/^hadron ~ 2. Instead, as one can 
easily recognize from the results presented in Figures 7.15 and 7.16, the multiplicity 
ratio takes values in two different regimes. In the "band" of values 1 — 5 for the case 
of the lower first impact and 100 — 160 for the higher impact. The larger values 
of the band in this latter case are justified by the fact that the shower is far more 
developed, given the altitude of the impact, and therefore is characterized by an 
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even more dominant electromagnetic component. The energy ratio, on the other 
hand can take small values, in agreement with the values observed in Centauros. 
However, notice that both the black hole and standard simulations show a complex 
pattern for these ratios and in addition they are characterized by large fluctuations 
for varying energy and number of extra dimensions. Furthermore, there does 
not seem to be a statistically significant difference in these observables between the 
benchmark event and the black hole mediated ones, at least for event heights greater 
than 500 m above the detector. We conclude that (a) either these observables may 
not be suitable, especially given the limited statistics of the existing and future 
experiments, to discriminate between black hole mediated events versus standard 
ones, or (b) that the differences in these ratios appear in events with initial impact 
in the range — 500 m from the detector. 

• Multiplicities 

For the multiplicities themselves the situation is much cleaner. In Figures 7.17-7.22 
we show the behaviour of the total as well as of some partial multiplicities in black 
hole mediated events and in standard events as a function of the energy and for a 
varying number of extra dimensions. 

The curves are very well fitted in a log-log plot by a linear relation of the form 

N = 10 9 < n >f£ (n) (7.87) 

with intercepts q(n) and slopes cr(n), that increase with the number of extra dimen- 
sions n. We present in Figures (a) results of the simulations performed with a first 
impact taken at 5,500 m, while Figures (b) refer to a first collision at 15,000 m. In 
Figures of type (a) the slopes of the benchmark events are smaller than those of 
the black hole events and show a larger intercept. This feature is common to all 
the sub-components of the air showers. A simple explanation of this fact is that at 
lower value of the impact energy, the number of states available for the decay of the 
black hole is smaller than the number of partonic degrees of freedom available in 
a proton-proton collision. We recall, as we have already discussed in the previous 
sections, that our benchmark results define in this case an upper bound for the total 
multiplicities expected in a neutrino-proton collision. Therefore, in a more realistic 
comparison, we would discover that the black hole and the standard results should 
differ more noticeably. The large multiplicity of the states available for the decay 
of the black hole dominates over that of a standard hadronic interaction, and this 
justifies the larger multiplicities produced at detector level. As we increase the al- 
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Figure 7.15: Ratio between N em (total multiplicity of photons and e^) and N had (total 
multiplicity of everything else) as a function of E\. The first interaction is kept fixed 
at 5500m (517 g/cm 2 ) (a), or at 15000m (124 g/cm 2 ) (b), and the observation level is 
5000m (553 g/cm 2 ). We show in the same plot the benchmark (where the primary is a 
proton) and mini black holes with different numbers of extra-dimensions n. 
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Figure 7.16: As in Figure 7.15, but this time we show the ratio between E em (total energy 
of photons and e^) and E had (total energy of everything else) as a function of E 1 . 
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Figure 7.17: Plot of the total particle multiplicity as a function of E ± . Case (a) is for an 
impact point of 5,500 m and case (b) for 15,000 m. 
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Figure 7.18: Plot of the multiplicity of photons as a function of Ei, (a) 5,500 m, (b) 
15,000 m 
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Figure 7.19: Plot of the multiplicity of e as a function of E±, (a) 5,500 m, (b) 15,000 m 
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Figure 7.20: Plot of the multiplicity of \i as a function of E±, (a) 5,500 m, (b) 15,000 m 
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Figure 7.21: Plot of the multiplicity of n~ as a function of E±, (a) 5,500 m, (b) 15,000 m. 
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Figure 7.22: Plot of the multiplicity of protons as a function of E\, (a) 5,500 m, (b) 15,000 
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titude of the impact, in plots of type (b) we find a similar trend but the differences 
in the total and partial multiplicities are much harder to discern for black holes and 
benchmark events. In fact, for collisions starting at higher altitudes the showers are 
all fully developed and the differences between the two underlying events are less 
pronounced. 

Another feature of the black hole events is that the slopes and the intercepts of the 
various plots, for a given choice of altitude of the impact, are linearly correlated. 
To illustrate this point we refer to Figures 7.23-7.24 from which this behaviour is 
immediately evident. To generate each of these figures we have plotted the param- 
eters (a, q) of a corresponding plot - for the total or for the partial multiplicities 
- independently of the specific number of extra dimensions. The results shown in 
these figures clearly indicate that the relation between the intercept q and the slope 
a appearing in Eq. (7.87) is linear and independent of n 

q = aa + (3 (7.88) 

with a and f3 typical of a given setup (photons, total multiplicities, etc.) but 
insensitive to the parameter n. Therefore, black hole events are characterized by 
particle multiplicities on the ground of the form 

N ground = 10 a ^El (7.89) 



• Lateral distributions 

In Figures 7.25-7.29 we illustrate the results of our study of the lateral distributions 
for the total inclusive shower and the various sub-components as a function of the 
incoming energy E\. The average opening of the shower as measured at detector 
level is plotted versus energy in a log-log scale. Notice a growing opening of the 
shower as we raise the energy of the black hole resonance, which is more remarked 
for a lower number of extra dimensions. In contrast, the benchmark simulation 
shows a small decrease (negative slope) with energy. The larger opening of the 
shower in black hole mediated events - compared to standard air showers - is due 
to the s-wave emission typical of a black hole decay, which is very different from an 
ordinary collision. Contrary to the case of multiplicities, here simulations of type 
(a) and (b) show a similar trend, with very distinct features between standard and 
black hole events. Notice that in this case the difference in the partonic content 
of the two different events (benchmark versus black hole mediated) is less relevant, 
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Figure 7.23: Parameter fit for the intercepts and the slopes of the curves in Fig. 7.17 
for the total multiplicities. The numbers over each point in this plot indicate the value 
of the extra-dimensions. The (a, q) parameters are fitted to a straight line q = a a + (3 
independently of the numbers of extra dimensions, (a) is the fit for 5,500 m, (b) for 15,000 
m. The benchmark is also shown in the plot, but has not been used in the fit. 
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Figure 7.24: Paramater fit for the intercepts and the slopes for the curves in Fig. 7.18, 
now for the multiplicity of photons, (a) is the fit for 5,500 m, (b) for 15,000 m. 
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since it is the geometrical fireball emission in the black hole case which is responsible 
for the generation of larger lateral distributions. 

Also in this case we discover a linear relation between average radius R of the conical 
openings and energy, relation that can be fitted to a simple power law 



In analogy to Figures 7.23 and 7.24, we show in Figure 7.30 that for a given setup 
there is a linear relation between slopes and intercepts of Eq. (7.90) 



with a' and f3' typical for a given setup, but again independent of n. 

7.8 Discussion 

A rather detailed analysis was presented of some of the main observables which charac- 
terize the air showers formed, when a high energy collision in the atmosphere leads to the 
formation of a mini black hole. We have decided to focus our attention on the particle 
multiplicities of these events, on the geometrical opening of the showers produced, and on 
the ratio of their electromagnetic to hadronic components, as functions of the entire ultra 
high energy spectrum of the incoming primary source. We have shown that in a double 
logarithmic scale the energy vs multiplicity as well as the energy vs shower-size plots are 
linear, characterized by slopes which depend on the number of extra dimensions. We have 
compared these predictions with standard (benchmark) simulations and corrected for the 
energy which escaped in the bulk, or emitted by the black holes at stages prior to the 
Schwarzschild phase. Black hole events are characterized by faster growing multiplicities 
for impacts taking place close to the detector; impacts at higher altitudes share a simi- 
lar trend, but less pronounced. The multiplicities from the black hole are larger in the 
lower part of the energy range, while they become bigger for higher energies. We should 
also mention that, given the choice made for our benchmark simulations, here we have 
been considering the worst scenario: in a simulation with an impacting neutrino it should 
be possible to discern between the two underlying events, whether they are standard or 
black hole mediated. The lateral distributions appear to be the most striking signature 
of a black hole event. Due to the higher prs involved, they are much larger than in the 
benchmark standard simulations. 

Our analysis can be easily generalized to more complex geometrical situations, where 



R = W'^E^ {n) . 



(7.90) 



(7.91) 
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Figure 7.25: Plot of the average radius R of the core of the shower of photons as a function 
of E\ for a black hole with a varying number of extra dimensions. The benchmark result 
is also shown for comparison, (a) 5,500 m, (b) 15,000 m. 
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Figure 7.26: Plot of the average size of the core R of the shower of e as a function of 
E x . (a) 5,500 m, (b) 15,000 m. 
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Figure 7.27: Plot of the average size R of the core of the shower of \i as a function of 
E u (a) 5,500 m, (b) 15,000 m. 
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Figure 7.28: As above for tt 
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Figure 7.29: The proton core size as a function of Ei, (a) 5,500 m, (b) 15,000 m. 



234 



7.8. Discussion 





Figure 7.30: Parameter fit for the curves in Fig. 7.25, describing the openings of the 
showers of photons, (a) is the fit for 5,500 m, (b) for 15,000 m. 
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a slanted entry of the original primary can be envisioned and, in particular, to the case 
of near horizontal air showers, which are relevant for the detection of neutrino induced 
showers. These are characterized by a larger cross section compared to the vertical ones 
and therefore are more likely to occur. The strategy presented in this analysis, which 
is limited to primaries entering vertically, does not change substantially in this more 
general case, except for the geometry which should give an asymmetric opening due to 
the directionality of the event. However, the main physical properties of the showers 
with an intermediate black hole should remain unchanged. These characteristics, in fact, 
are not sensitive to the geometry of the direction of the event but are due to the larger 
multiplicities of the decay of the black hole resonance that forms after the impact and to 
the s-wave structure of its instantaneous decay. 
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